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ABSTRACT 


This thesis examines stochastic linear programs, 
i.e., linear programs in which some of the elements 
Ofte ;cOcliicCrentymatrix, resource vector, Or price 
vector are random variables. The three classical 
approaches to the stochastic linear program, linear 
programming under uncertainty, chance-constrained 
programming, and stochastic linear programming, are 
reviewed. The emphasis is on stochastic linear 
programming and the simulation of stochastic linear 
programs. A new estimator is proposed for the decision 


vector in a system of stochastic linear equations. 
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CHAPTER ol 


AN INTRODUCTION TO STOCHASTIC LINEAR PROGRAMS 


Aystochastic linear program is a linear program 
in which some of the constants are random variables. 


Thus if a stochastic linear program is represented as 


Zo eA Se Cee X (Ea theal)) 
= 
subject to 
IN See 18 (dig les) 
and 
sw F Ci l—3)) 


Whevew Am edie OV elie ad tar manda, eC, aNd ix sdreuvecrors 
in m-, n-, and n-dimensional real space (Rew. Rei and 
Re"), respectively, then some of the elements of A, b, 
aie ecetiave  cenancovmvaria)| eS mem neaestochastic: lunear 
program we must also assume that the joint distribution 
OLA; meena Cais knOwn,awheseas sin -asclhassical lanear 
program we assume that A, b, and c are known constants. 
Naslund (1965) is the only survey on stochastic linear 
programs, and while it provides a good introduction to 
Eneesubjectmancd acOntains, a fairly extensive bibliography, 
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will provide the missing details of the techniques 
involved and illustrate the techniques with simple 
examples. 

iiex is cheoceswbefore*the valdestof Ay \b; and 
CROCGoner nOWnymcoensit Us possible tehat “xomay ‘not 
Satighva) -l-2) @atebomek cand cb become known; Le dpa 
may be infeasible. There exist two methods which 
deal with this infeasibility problem - linear program- 
ming under uncertainty and chance-constrained program- 
Ming. These two methods are discussed in Chapters 2 
and 3. Chapter met includes the details of using 
the decomposition algorithm in linear programming under 
uncertainty which was first suggested by Dantzig and 
Madansky (1960) but has not been elaborated elsewhere. 

li taeeewchosen after*the values of the random 
WVeetaDleés@epemobserved, then svand z* are also random 
Svariables @melne study of the distributions of x* and * 
is referred to as stochastic linear programming. 
Stochastic linear programming is based on the solution 
of a system of stochastic linear equations, just as 
linear programming is based on Mieuaomer ken of linear 
SequationsmeeaA first order estimator is proposed jin 
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alternative to an estimator introduced by Tintner ((1955) 
p. 204). In Chapter 5 we will discuss stochastic linear 
programming and offer a suggestion for extending these 
estimators to stochastic linear programs. Further, since 
Monte Carlo simulations have been suggested throughout 
the literature as a method of analyzing stochastic linear 
programs, Monte Carlo simulations were used to analyse 
the examples in Chapters 4 and 5. As far as the author 
is aware, the majority of the material in Chapters 4 and 
BeLseorrgiiial.. 

The emphasis in this thesis is on the simulation 
GOrestOctastics lilear Cguations and programs. It was 
originally hoped that these simulations could be done 
in APL. However, the long response time and space 
limitations proved to be too great a handicap. There- 
fore, after testing the APL random number generator in 
APL (see Appendix 1), it was translated into an assembly 
language subroutine (see Lewis, Goodman, and Miller 
(1969)) and then was used in conjunction with a FORTRAN 
program to simulate the stochastic programs (see Appendix 
2). We also attempted to use APL notation as a substitute 
for the conventional mathematical notation in Appendix 1 
Dy writing all the formulae in both mathematical and 
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CHAPTER 2 
LINEAR PROGRAMMING UNDER UNCERTAINTY 


Zeiaeerneroaquccion 


Linear programming has proven to be extremely useful 
Ineangreat variety of applications. However, the formula- 
tion of a linear programming model tends to be very 
Testrictivesandserigid.paOne reason ,for.this is the.deter- 
Ministic nature of the classical general linear programming 
model. Usually a linear programming model may be written 
assLollows: 


Find 


n 
max are Cee, (C2r eh) 


n 
EPO a er ee 15.2) ea. pt (2.1-2) 


and 


ia SSO) Te Flare ipretet ens Live (221-3) 


Max, 6c 9X (2.19) 


subject to 
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5 Seg ; (21 =3 4) 


WoCLemCe omd0e ne OVE lEVCCLOM  sOe1 Seanemebyel avector, gAliS 
aneiebyvanemtatrixsandese1sS an n by levector. «The model 
Cer redmoyven er — te clas) anda — saw | be. tTe— 
ferred to as the deterministic linear programming model. 
The assumption that c, b, and A are known constants 
1S notlsalwavs justziied. In practical situations the 
values of these constants may not be known and they may 
be better represented as random variables. This lack of 
complete information is not enough in itself to make the 
deterministic model any less valid, as it still may be 
Deere Omi oOlvcmUacmUunGCEtEdiNtelos DV alaking ssOmeuad @OC: 
SpeChitCats Of -mOunC Deals Agand athnCne COmsOlVveslLhcacdaver— 
Mvovsttecemodel a s1iis method of replacing, the stochastic 
elements by constants may or may not be justified depending 
CONeLUemScChStGLyicy eOGbeLies Model to, Vabmatlons in the, para— 
meters. If this model is insensitive to parameter changes, 
then a deterministic approach may be justified, since the 
introduction of uncertainty into the linear model makes 
thesanalysisymuch more complicated. Stochastic models are 
inherently more complicated than the deterministic ones » 
Pome BUG s ee LoCeLedcONS ween oouly ys liethe .slochastic 


Modem use COstLake Into wccounteall jnossible variations of 


f*E=1.8). —, £2 
v 2 A .todesv fo yt m m6 -aiod jxotosy I vd & os 42 So 
fehow at .zotosy f yd h ns @t x bn6 xlasgm A ye at th ho 
oes ed ifiw ("E+2:5) bas ,('S-5.f) , (*I-t.e) yo Sock 

S : 7 ~ 

. lebow Onigiens el 3e5Gi I Ss terormsteh sis ee og boxe}, : 


a 
SsKe2en0S awodd ate A bis do JE07 Koktqmada set | oes 


\ ‘@n3 eeoldsentice [e0ita6i9 ol -beltiteot ayewia ton ei 
‘Yam ¥oss fins twoad ed ter yon aines@noo sgecs te eealay - 


So Stoel ein? ..asitéiwsy mobdar es botnsaetaey esse i 


of3 stem OF BMeedi ni dovons ton ei neldamricial sxielqms! 


ea 


od yom ticve 37i 26. .filav sacl vus febonr si sega tessen 


* . 


ned Bs smpe perivianm yd esitiisdcrocny wit oviceet- of sesaeeee 
——— - — a = 
ba — 


~ Taieb en? avibe oF feds brs A Ore a ‘2 20 ants aes a 

pisastivots e639 grisel[as: to Borcem aid? ~.lethem viselnin a 
mifbaegen. Beizizedt sd feb yan th vem ear edetos vd atosmale,) 
-s269 Git ci wnolyaitay as [obom site Yo: ytivisctense edgy 0 | 


@epned> se7etmisg ct avitiensuai =i_ lobam eid 2 


eaves isGom ragni] sc iit vieiseracke to apis 
Sth eienom. Jitesroess betes} fonss son Ain ninetehale 
Geno Sitdelnivwesed sds wes bsvaablgmar Sion wiht 
aitasioese ots 22 .yfoexta -2enensx - aol 9 
40 anobse) zie Bidteseq ifs Jayenge arm ie 


7 - : : 


f ~ 
on 7 : 7 : _ — 
. Na = _- _ au 
; ' 7 im ao 7 - 
: > ao _ to 7 
t ._ = a -_ a a 


the problem, then the equivalent extensive form of the 
deterministic model will be large. Secondly, the objec- 
CiVvemeunctmlOongsoOr vdestochasticumodely 1s stbypicalis., non 
LINeaAg. MLAS y,, sincesthe distribution, off cacharandom 
element must be specified, more information is necessary 
to specify a stochastic model. Since a stochastic model 
may be more realistic but more complicated than a deter- 
mMinistic model, it is necessary to make a trade off between 
reality and complexity. 

Dantzig (1955), and Ferguson and Dantzig (195¢) were 
Lien LSoiLOmCOneotders Stochastic, linear programming models, 
andwabpliece their modeleto thesgallocation Ofaaixcramts to 
different routes subject to uncertain passenger demand. 
Their model may be stated as follows: 

BindmasvectOlexa to Mingmize 

min E(e'x ae EM y | (2.1-4) 
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ditienslonal vector yis chosen to™solve the»equation (2,1-5). 
Dantzig (1955) also gives a computational procedure for this 
two-stage linear programming problem. The probability 
denver lymrunct ron Of Dewill be=assumeds to pe=known. This 
problem could be referred to as linear programming under 
risk since the probability density functions are always 
known. However, the more usual name of linear programming 
under uncertainty will be used. In practice the derivation 
of the appropriate probability density may not be obvious 
and may have to be obtained from relative frequency data 

oer) subjective probabilities. 

Dantzig and Madansky (1960) use the dual of the 
linear programming model under uncertainty to obtain a trans- 
formed model to which the decomposition algorithm may be 
applied. This paper and Dantzig (1955), referred to 
previously, are the main sources for the material presented 
in the remainder of this chapter. 

Madansky (1960) has demonstrated that the substitution 
Ofsthe expected values for the random variables does not 
necessarily provide an optimal solution. However, such 
a formulation can provide a bound for the optimal solution. 
This bound and others considered by Madansky are presented 
in Section 2.4. 

Managasarian and Rosen (1964) develop similar bounds 
to those of Madansky for the model containing stochastic, 


non-linear constraints. Wilson (1966) develops a generali- 
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zation of Dantzig and Madansky's method to obtain an 
approximate solution for the stochastic, non-linear 
model. These techniques involve non-linear programming 


Metnoasmend aremnou conslderea ine this thesis. 


2.2 Linear Certainty—-Equivalence Theorem 


Consider the following formulation of the deter- 


mMinistic, linear-programming model: 
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We will assume that the e248 are random variables with 
known distributions, and that the Sas must be chosen 
before the aes aGemoUSeLlyed. ASeaeresULi (2 .2—1)s Cannot 
be optimized since z is a random variable. However, we 
May replace z by its expectation 
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the expected value representation, it will be used 
exclusively throughout this thesis. 

To illustrate the above type of problem consider 
the following example from Dantzig (1955). Suppose a 
dietitian wishes to formulate a minimum-cost diet. 
Since the prices may be subject to fluctuations due 
to changing supplies and demands, they are not known 
with certainty, and the deterministic model is not 
applicable. We shall assume that prices are random 
variables with known distributions. .The dietitian must 
prepare a diet which not only meets specified nutritional 
requirements, but also has a minimum expected cost. Let 
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hp be’ the number of units of the 4° food purchased at 


a Unit ‘pricervoL £0 ae be the number of units of the 

sa Niearents contained Tn asunut “Ofe the Be food, and be 
be the minimum requirement for the aes nutrient that is 
Fequnmeds Dye almindivrdual®.- BThewproblem cane then? be 
formulated as selecting the aS to minimize the total 
expec eechicese™ (22-4) sisubsiec teeo tunes nutri onalarequire— 
iMeiTiuce (ino —2) Band Li2.2—3) 8 Belhise problem cane be-solved 


by applying the Linear Certainty-Equivalence Theorem 


(Wagner (1969)): 


Theorem 1 : Linear Certainty-Equivalence Theorem 
‘Assume that ais and S ica orm tipng (2.2 — 2.) 


and (2.2-3) are known, and oe are random variables 
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* 
prior to knowing the exact values of a then ea is 
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SUubpecomeOm co a mrandan(2.2-5)). 

To solve the stochastic diet problem, the expected 
values of the prices need only be substituted for the 
prices and any linear programming algorithm applied to 
the resulting problem. Although the deterministic equi- 
valent form of the diet problem is simple and easily 
Oped ned, emsnouldgno tebe assumed) Chatetiis. 1S always 
SOnmee Ltt general, the expected value may not be substituted 
for the random variables and the deterministic equivalents 
to the stochastic problems are not usually as simple as 


the diet example. 


2.3 Two-Stage Linear Programming Under Uncertainty 


In the two-stage linear programming under uncertainty 
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E [min (Gilaax Choe fs! y) | ; (2.3-1) 


where c and £ are known nj7 and n-dimensional vectors, 
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respectively, and satisfy the constraints 
SS is Vs (223-2) 
xe OMe ¥Y o> 10 ; (25 3-3) 


A and B are known m by n, and m by n, matrices, and b is 
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an m-dimensional random vector with a known distribution 
function independent of x. In the first stage we select 
tieevalucsvom xeknowmgsonly the Wwalues of A, B, c,) and 
Peo LCT moc LeCLi nga x, swe Observe —aurancdom Vector b.. In 
the second stage we choose Vise Opsa cist yea 2 )a-pagt he 
vectors x and y are selected to minimize the objective 
funebion «ln general, the, solution will have a, value’ of 
See Clie Wilsl aD GeWoetser Coat lescrOre thawb jobsenved but ta 
complete strategy or contingency plan is needed for y, i.e., 
for every possible b observed there is a particular y which 
will be used in order to minimize the objective function. 
iiasaomeapDatLent, that, the solution to ethis type of problem 
is more complicated than the deterministic one. 

In. the above formulation it is necessary to add the 
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Ofebye4 -c- pe giventanyiteasible x and a possible b, the 


set of linear equations 


By = b - Ax (2.3-4) 


must always have a non-negative solution. Equivalently, 
we may assume that there exists a convex set K such that 
iexse Koeethen x) 270" and there exists a y for each pb 
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subject to the constraints 
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The vectors nN, are unrestricted tn sign. “Tt 
Ny ret eye OF (255-1 1) 


then the problem becomes 


L ] 
rer (Py bi ¥, + Po bo Vee ae Any Sp Py he ME gy He eked) 
oral 


subject to 


§ i i] 
Py She YoetePoe® @hoiterehitePyiA’ Yas ¢ 
i] 
ees < 
B' Y, <—£ (2.3513) 
@ . 
@ e 
\) . 
Be YW< i - 


Except for the being unrestricted in sign, the problem is 
now in the form required for the application of the decom- 
position algorithm. The first constraint represents the 
major problem, while the last N constraints can be used as 
a Single subproblem as Dantzig and Madansky (1960) 
suggested. 
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can be considered as a subproblem. This approach may be 
extended readily to include problems where B is random. 
165i B. is the random matrix corresponding to bis then the 


aia subproblem is 


i te (2.3-15) 


The following analysis is more general than that given by 
Dantzig and Madansky (1960). 

In the original dual problem, the dual variable Nn, 
are unrestricted. In order to apply the decomposition 
algorithm or the simplex algorithm the variables must be 
non-negative. If we let 1, =a, - B,, where a; 2 0 and 
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simplify this problem let us define the follow- 


and matrices: 
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is necessary to assume that each subproblem is a bounded 
set; if not, certain modifications can be made to enable 
aesOLucLOneeco Her round (Dantzig and sWolte(1960))". = The 

decomposition algorithm can now be applied to this model. 
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Bit <“p. t=“and 1, 2 OF, tater. ey. (23 0) 


Moreover, the representation of 1; in terms of the extreme 
points of S5 can be substituted into the model (232 3)8, 
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we must now solve for K variables Oj. and 1;) where 
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If the revised simplex method is used, it will not be necessary 
to find all the vectors: in each W; Since there are only 


n, + N basic vectors in the basis at anyone time. 
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(Cp) g = b. Tee (2°.3=49) 
and 
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Once the transformation is complete, the optimality criterion 
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and 


to wae 2s, ENE TISt. be solved. 

It is obvious that a great deal of work is required 
to solve such a stochastic model and that for small problems 
Leglce probabil yaletterm to solve =the smodel 9(2.3-6), (2-3-7), 
anom 220-0} sand ereserve the*dualmtormulation for large 
problems. In any event it should be clear by now that 
stochastic models should not be used indiscriminately and 
DEC eVonecr pews nelpeto clarity the particular class of 


problems in which they should be used. 


2.4 Example of Two-Stage Linear Programming 
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This problem could be represented as 
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subject to 
PO Xe te By = Daa 
or equivalently 
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variables Y} and Y5 are, respectively, the number of ship- 


MenUCSMCO elim croMmetcneswarehnouSe ateaLunit cost £. and £ 


i Ze 


respectively. 
tiewevacsignecome antegral values to c, £, A, and B, 
we obtain the following sample problem: 


Z = min E [4x, + 19x 
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Suppose now that b can have one of two values by = mOUL, 45)' 
or b, = (36, 54)', with probability P, = 0-25 and p,= 0.75 
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TADDERee =) 


Second stage strateqy for y 


Number of Shipments 
Demand 
aie 


(30,45) 


(367.94) 


thus) Stratedy aamplies: thatyit by occurs then there is a 
surplus which is shipped to the warehouse and the resul- 
PAntmCOS us 9220.85 Liethe demand bo occurs then there is 
a shortage, and the factory must bring in additional 


supplies from the warehouse and the resultant cost is 
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2.5 Inequalities for Linear Programming Under Uncertainty 


In this section we shall discuss a few alternative 
models ye cie linear programming problem under uncertainty. 
These alternative models will cen be used to show the 
influence of uncertainty in linear programs. Madansky (1960) 
was the first to use this type of analysis and a good summary 
of his work can be found in Wagner (1969). 
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must be made before all of the information about the system 
is known. For convenience, let us denote the objective 


function by 


min E [C(b, x)].= E[here and now] (2). 5-1) 
x 
where 
(ONS, eS Et ule Me oa se tee (225-2) 
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The following analysis depends on the assumption that for 
anya lve ex wd ncmOmLlloresalwayceexists a feasible solution 
yy to the here-and-now model. The distribution of the random 
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EORaD IO Lys. | 
iivaggdecision cn) the variables x, and y is postponed 
until the random elements are observed then the model is 
called a "wait-and-see" model. The problem then is to find 
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pUEE 
SUbjec mt Omesno Oy medniOs.2. 5-4) wheresh) 1s Known. 9\lhesexpected 
value of the optimum solution for this model is 
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optimum value, weight these values with the corresponding 
probabilities, and then calculate the expected value. This 


would involve solving deterministic linear programs of the 


form 
Tere (Ca Day) (2.5-5) 
Rin 
subject to 
NBS ap WP be (2-5-6) 
and 
ns SMe p i Ow 7 (26 eu) 


Wht Cut ler2 a heres aN 
Consider, for example, the problem in Section 2.4. 


If all possible problems were formed and solved, then 


E[wait and see] = ($180) + 2 ($216) 


$207 
Note that E[here and now] = $224.50 and hence 
E[wait and see] < E[here and now] . (2.5-8) 


Madanskye (900) mhas shown that (2-5-8) 1S generally true. 
Since the wait-and-see model is not subject to uncertainties, 
LUeLomoesyocem OL COnplActe,intormation, and thus it 1s reason-— 


able that the expected minimum total cost is less. 
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If the average values of the random elements are 
used in place of the random elements then the model is 


called the "average-value" model and may be stated as 


follows: 
Tn Cue Xeeoecuey eo [average ava lucl| (2.5-9) 
lied 
subject to 
EN pete 1) ye 1h. [isi (255—1,0) 
and 
Kies Ol, Vara Oe ee (22) 


If we consider the example of the previous section, 
thengthne optimum value for the average-value problem is 
E[average value] = $207, which is equal to the solution 
to the wait-and-see model. Madansky (1960) found that in 


general 


E [average value] < E[wait and see | (2.5-12) 


andsinvpaneiculan estates that: 
PieCliCm DL Ola) Milt Ve OmmeaslinewOnmticsselponsb'S 15 
COUNtaALbLYsaddalive or the set of DeseicewiLthneproba— 


bility one finite in number, equality holds between 
E[wait and see] = E[average value] (2.5-13) 
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Ifthe optimal “feasible value of x in the. average- 
Val Uemmooc Peo axa ne thengthe model obtained by using x* 
acetieerirse—stage strategy and y as the second-stage 


strategy is called the "first-stage, average-value approxi- 


mation". Mathematically this means that 
Mine (Ciel by )= = Scere rami n (ily) (225514) 
ey | q 
subject to 
Bays = Die eae ne (225 —5)) 
and 
y > 0 s (2755=16;) 
Hence, 


Ec xe + asa ey) = E[first-stage, average-value approx. | 


a (25 


Again the expected value of this approximation can be 
FOUNnGmSCINO Ly eeOGearG1SCretC distribution: Lunctron by solving 
all possible problems as was done in the wait-and-see model. 
In the example, the first problem would be to find 
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and 


Clearly, there is no feasible solution to this problem 
Since it violates the assumption that there always exists 
a y which is feasible for each pair (b, x). 

It may be shown that, the expected value of the 
first-stage, average-value approximation, if it exists, 


is an upper bound for the expected value of the here-and- 


now method, i.e., 


E[here and now] < E[first-stage, average-value approx.] . 


(2eS—15) 


There is an interesting economic interpretation to 


the difference 


E[here and now] - E[wait and see] 


2) 
= 2 Ot Om 20) B= 51 50 


Iimecndielt ere eleceEsathe COS@ Of UNCertainty. This is the 

GxELagcostediuc toncnoosing), X prior to observing the random 
demandeope ut isealso the amount the firm would be willing 
Loma MEOmLeCmol omeOnDLecaleL the demand perfectly. Suppose 


that the difference 


E[first-stage, average-value approx. = E[average value] 


(225-20) 
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exists where both of expectations are simple calculations. 
If this difference is "small" then it may be more economi- 
Calelouuse x*sas thnesfirst-stage optimal decision variable 
for the here-and-now model. This plan will be more 
economical if the extensive form of the here-and-now model 
is very large. Thus the difference (2.5-20) indicates the 
closeness of the optimal values for the average-value and 


the here-and-now models. 
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CHAPTER! 5 


CHANCE-CONSTRAINED LINEAR PROGRAMMING 


Sewer Ormulation Oo: the: Model 


In this chapter we shall consider an alternative 
model for the classical linear programming problem in 
whicheithe coefficients are, subject to stochastic varia- 
tions. It was developed originally by Charnes, Cooper, 
and Symonds (1958) in considering the scheduling of 
heating oil production, the demand for which is subject 
to marketing and storage constraints. This model, some- 


times called the E-model, may be stated in simplified 


form as 
ine ECs x | (3.1-1) 
x 
subject to 
Pra (Aa xem) Oe, (so?) 
fee 
n 
ihe le) au a > b;) > Oss P1062, eee Ny 
j=1 
where 
OR< ae Sak ' i=1,2, /M , 
and 


a WY 2 (s24—3) 
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The deterministic constraints of the classical model are 
replaced by the chance constraints (3.1-2). In the 
solution of this model, which will be represented in 
Section 3.2, the chance constraints are replaced by 
equivalent deterministic constraints using the dustri—s 
Bier one UNer Ono tao me coarnesmandeCocperg( |I59)e give 

a more general and theoretical treatment of the above 
type OLeproblem,. 

The problem of scheduling SneOLUetankeiaesrleetnis 
considered by Charnes and Cooper (1962). An optimal 
Ship-chamecringme0Old cysis  LOUNGeSUbyeCL tO ane uncertain 
shipping demand. The cost of chartering ships is assumed 
to be a random variate which depends on a normally distri- 
buted shipping demand. As in the heating oil problem, 
they are interested in a scheduling problem, consequently 
they introduce a linear decision rule into the model (see 
SECClOnso a5) )a 

Charnes and Cooper (1963 a) consider, as well as the 
E-model, some alternative forms of the objective function 
rl Pee eune Objeculverrunction 2S replaced by 


Hina VOM Ke My 7 Cc. x val j (3.1-4) 
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the V-model. The model formed by replacing the objective 
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is called the P-model. 
Miller and Wagner (1965), and Symonds (1967, 1968) 
considerereplacingethesym constraints) (341-2) ) by-ai single 


JO2 DeESCOonstlLaint 
Prem (Axes) Db). > 0 : (3.1-6) 


This model is simpler since only one value of qa has to 
be specified, but the solution is more complicated. Let 
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with the sign of ny unrestricted. The function Fm) 

is nonlinear and, according to Wagner (1969), is Se Ay 
convex. Thus nonlinear programming can rarely be used 
as the problem in its present form. However, if we take 


thesnavuuna lee ogamiuehm= of eboth ssades: 10 fae3% L-10)i. ieus, 


lms 


ln F. (nj) be Eniscl a, (@3.e—] 1) 


then the left-hand side is a convex function when be is 
a normal, gamma, or uniform variate. Consequently, the 
joint chance-constrained program can usually be solved 


using nonlinear techniques. 


3.2 Solution of the Chance-Constrained Model 


The chance-constrained linear programming model may 
be used in place of the classical linear programming model 
when it is desirable but not absolutely essential that the 
GonstGaimts hold. = lhe 05 in (3. b—-2) Sarre  speciivedsconstants 
in the interval [0,1] and are usually chosen close to one. 
The quantity l-a; represents the allowable risk that x5 
will not satisfy (3.1-2) once the values of the random 


variables are observed, i.e., 
n 
i” Gel egy eae oe UE . (3.2=1)) 


Thus it is necessary in the chance-constrained program to 


define the vector x as feasible if and only if x satisfies 
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(Srl ve ctiGms Seto re Tne. Obj ectinvve is to select the 
‘beste nonnegative solution that will “probably” satisfy 
EiescConSccCramncspn xe Darter whe values Of the random 
Variables A, ~b and Gc are observed. 

Let us consider a particular subset of the class 
of chance-constrained programs where as. and a; are 
FknOWne Constants eine Variables ec. and Oe are random 
and statistically independent for all i andj. Again 
ChesAUSt ri DutToOnerUnee Ton, F, (b) = Pr (b; <a) eLOm cach 
bs is known and is assumed to be normal with mean Hy and 
standard deviation O.- TOSsOlve=thise specific chance— 
constrained linear program, the probability relations 
must be replaced by equivalent deterministic ponee cen nkss 
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is equivalent to 


(332-5) 
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However, the right-hand side,of this last inequality would 
be more useful if it were expressed in terms of a standar- 


Qizved' normalavariate, i.e., 


-1 _-l 
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where 
on x 
wat] 719 
Eves 6(z.) =f — e abe (Cee) 
—0o V2T 


th 


Therefore the i chance-constraint in (3.1-2) can be 


replaced by, 
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and the chance-constrained program defined previously in 


Section 3.1 is equivalent to the classical program 
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min [4x, + 12x, + 9x5] 
xs 
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than the two-stage linear program under uncertainty. 
Chance-constrained programming only indirectly evaluates 
the economic consequences of violating a constraint, and 
it makes no distinction between a small and a large 
shortage of production. Also the values Os of the 
chance-constrained model must be specified before a 
solution can be found. 

However, the chance-constrained formulation does 
have some advantages. The equivalent classical program 
for the chance-constrained model has the same number of 
constraints as the original program. Whereas the greatest 
handicap of the two-stage linear program under uncertainty 
CP OmeCtoc uc Como) McmClCe Les intal Ces i Zest lt Smequivalent 
Classica lm pmOog baller UGtietsthie SOlLuULLOnMOm the chance— 
constrained model is much simpler than the complete 
strategy solution for a linear program under uncertainty. 
Finally, the chance-constrained model requires less in- 
LOLMationsalotltmenc distribution of b.- only the as 
fractile: of be need be known. Therefore these advantages 
and disadvantages must be weighed carefully when deciding 


on a particular model for a given problem. 
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CHAPTER 4 


SYSTEMS OF STOCHASTIC LINEAR EQUATIONS 
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A stochastic linear programming model has the 
form 


Ze=si Maxey Cl xX (4.1-1) 


Ax <b (4.1-2) 


and 
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where A is anm by n matrix of random variables, b is 
Anesi._C Mens tenauE.eehOreotearandonayeriables mandgcyis 
an n-dimensional vector of random variables. The joint 


probability densmeyarunction of A, b, and’c, 
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the stochastic programming may be used to select the 
ebest" probability density function, defined to be that 
density function which maximizes with respect to 9 some 
given preference functional, e.g., 


meaxee te | 2m or | (4.1-6) 
GE 


Instead of the expectation, the mode, median, variance, 
skewness, or some other statistic of the distribution 
Opmzmiay sOesuSccmI Tm 41-6) an Ldeally, the solution, 
i.e., the density of z, should be derived analytically, 
buttastwell besseensthis 1s notepractical*for the 
majority of problems. Consequently, we must settle 
for some approximate solution. 

Suppose that the slack variables have already 
been added to the above model so that the constraints 
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and 


x =er 08D (4.1-9) 


Thessubvectors x4 and €, correspond to the particular 
submatrix Aj: Thus in the simplex algorithm we must 
ultimately solve the set of linear equations (4.1-8). 
Correspondingly, in the stochastic linear program we 


must solve a system of stochastic linear equations for 


the vdistri bution Obex, where 
AGSOR FOR: : (40-10) 


Pe CeO ee omc Ss SUMeC CMO uDema nel  Oyemematrieands b 
an m-dimensional vector both composed of independent 
random variables. It is the purpose of this chapter to 
investigate such systems of stochastic linear equations 
while the next chapter will discuss stochastic linear 
programming. 
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pendent normal variates, was by Babbar (1955). His 
apDEcoChmlsombasca one solving, Lor sxfusing sCramen s rule 
and then expanding the required determinants in terms 
GieticlmecG reactors. BHhe Chen ignores for no apparent 
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components of x as a ratio of two linear combinations 
of normal variates. (He claims that the ratio of two 
normal variates is approximately normally distributed 
if certain criteria are satisfied.) He then is able 
to approximate the distribution of this ratio by 
EranstOrnUngeelvsequdltitysintOordmVariateswhiirch wis 
approximately normal. It will be shown in the next 
section that this ratio of two normal variates can be 
solved analytically, hence no approximation is necessary. 
Babbar (1955) also extends this method to stochastic 
linear programming by assuming that only one basis may 
be optimal, i.e., the chance of some other basis being 
optimal is assumed to be zero. His approach to the 
problem has been criticized by Wagner (1958) as being 
restrictive and vacuous because he ignores the non- 


negativity requirements. 


4.2 Transformation of Random Variables 


There does not exist any literature pertaining 
directly to Phesanialyivcal solution som the stochastic 
Winecatmeduaclonssoteties forme(4 41-10). Most of the 
literature relates to systems of linear equations which 
abewsib jecuato stochastic errors, and to the determina-— 
tLOneOt error bounds for x. This type of error analysis 


is of no help when x is a random variable for which we 
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require a probability density function or moments. 

However, we can find an analytic solution to a system 
of stochastic linear equations by applying the theory 
of random variable transformations (see Hogg and Craig 
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The system of stochastic linear equations 
Ax = bos, (Ae) 


may be viewed as a transformation of random variables. 
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or equivalently 
xo b . (4.2-4) 


and select these additional variables so that the mapping 
is continuous, one to one, and onto then the inverse 
Mapping exists. Suppose it is possible to express the 


inverse mapping, b = A x ,.as 


2a Sid, = Sills) oi) Bea ero (4.2-5) 
and 


(x), k = mi + 4, and (4.2-6) 


If the first-order partial derivatives of the inverse 
functions’/are continuous, then we may express the 
Jacobian, 7, Or this transtormation as 
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This method while very indirect provides the only 
general method of achieving anpantal Visice @SOlUCLON ECO st he 
system of stochastic linear equations. In special cases 
the transformation approach may be shortened by applying 
theorems on linear combinations of random variables, and 
these particular cases are noted in Section 4.5. It 
should be obvious, that the transformation approach is 
impractical because of its complexity,especially for problems 
involving more than m+ m = 2 dimensions. 
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will have to be satisfied with estimates of the expected 


Value ;m=Vanianee, and Covariance of x" 
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The most widely used approximation for x is due to 
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To find the expected value of x the expected value 
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The only necessary condition for the estimator xn 
is that x must have continuous first and second order 
derivatives in the neighbourhood of the point 
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The two types of problems to which he has applied his 
approximation are either when A is random (Tintner (1955)) 
OpeDe se randomm uli nenereand Kaghavane('965))7 but note 
when both are random. 
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In the previous section we found the probability density 
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where A is an m by m matrix, and x and b are m-dimensional 
vectors. 
(1) 


The sequence of matrices A , A tou; wl Saesald 
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Next let us define the vector norm, ||x]|, of x 
as a measure having the properties: 
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Also we can use an at order approximation for (I-A) ~, 
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whereyIeassanim by m identity matrix. If the determinant 
ose Ja oy Mo 1s non-zero, then its inverse exists and 


we have 
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However, the probability of the matrix being singular 
Maye become non-zerosifythe determinant of the expected 
value of the matrix is close to zero and its variance is 
"large". Notice that |E [I +A Mo] |= 1 so we need 
only to keep its variance “small" for I +A Mo” to be 
non-singular. 


Westangnowesolvemror xfjin) (4.574) isinee weehave 


already assumed that Mo eXists; Laei7 


eS eee Me) be. (AMS 25) 
om (e) © 


If we can show that (A Mo*)” Osea Snnee orm nenmwercan 


app lyalieorcimomtOnig + A mot since the Pr{|I +A Mo* | 40} 


is one. From Theorem 4 we know that if ||(A Mo”) || <2 
vedio 4@) ai 


Siena g ee eee Meo Cele Lhe the Prt |) A we ly 


approaches one (see Section 4.7) then it is reasonable 


to apply Theorem 5. Thus by Theoremes (4. 5-90) sbecomes 


Me 7 


me Pere) eee?) 258 (425-6) 


x=M (ie eae 
—_ O (@) 
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se 


shandinndeh ajo 22. .xbzdim ht detent 1% ae “giteles 
. bee edeixe setevdi eft sods SISTA as eae + @ 


a, 


a ? ar 
(b-2, 8) - ao PEMA + ol) = eee 0 
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or equivalently 


ih 2 med 


b bat -1 
ee ee ieee (MNS A) Ge) i Mie Dee *(425-7) 


For convenience we will write (4.5-7) as 


GE SON eS pe te (4.5-8) 


a 


al 
fe) ° 


where M M 


From (4.5-8) it is now possible to define a first 


. “w~ 
OLOGr@eGStimadtOn Xerone~ as 


x= (I-MA)Mb=Mb-MAMb , (4.5-9) 


where the error ||x-x|| can be found by applying Theorem 
Cea NESmCoLIMacOnhmeusmialrsEeOrder in A ailevenein contains 
ChiewcLOsScmDNOCUC EAs eOaNex awe swill la LUnCmtne meaneand 
the variance fof "this first order “estimator. 

Applying the expected value operator to (4.5-9), 


we find that 


E [x] = © [mb] - & [MAM Db] 


Manan moe Mee |All Mer. (i| (4.5-10) 


feu He |b) bam andirecalD that’ & (A= s0.ee Rhus the 


expected value of x is 


fae | eee Ma (4.5-11) 


| ts ob 
(ane) od “Mt... - “+ A 


ae (T-2.0) stiaw Dfiw ow Somekisvinga aan lm 
ae i 

(842. b) dee. = My + a= 2) ee 4 
“i : we ai 
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= See 


S242) BS Gokteb ot Sliiescg wom ef 2k (6-205) meee 3 ‘s 
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Unfortunately, the-variance of x is not as easy 


torderive. First, by definition 


Exerc |Seoe ee: (4.5-12) 


: 


1X > 


var (x x') = E[x 


UST Nes 407) pet hee math 1x x x! becomes 


x x' = (Mb-MAMb)(Mb-MAMb)' (4.5-13) 
or 
x x' =Mbb'M'-MAMbb'M' - Mb b'M'A'M' 
+MAMbb'M'A'M' , (4.5-14) 


The expected value operator is applied to (4.5-14) to 


ODtal nwthiceet rst term Or the variance (455-12), 12e%, 
E[$ $']= ME[b b'] wm’ - M E[a] M Efb bt] 
-ME[b b'] M E[A'] M'+M E[AMb b'M'A'] M'. 
(425-15) 
Recalling that E [a] = 0, we can simplify (4.5-15) and 


obtain 


Bix x!) = ME[b bv Ma + ME [A Mb b* M'a'] M’ 
(a5 6) 


The expectation operator is an idempotent operator and 


thus (4.5-16) can be expressed as 
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Bex eee (Deo eA ME [beb'| MiaAt] M! 
(455-7) 
We) could use (4.5-17) and (4.5-11) to express the variance 
(Amb We butetnitsswoulrd be Of fittle- use. awe need to 
express the Var (X X') as a function of the known means 
and variances of A and b in order to calculate the value 


of Var (xX re To do this the following lemma is necessary. 


Lemma: Let A = [aj 5] be an m by m matrix of independent 


random variables such that 


i = © (4,5-18) 
and 
fav = 6 = [o$5] Bh, Soa) 
If M = [m5 5] Vopaneieyallemd bri OlecOns tants, ms then 
Ena Men | @=—en (Ss ee (4.5-20) 
/ 


& is a diagonalization operator which transforms the m- 


dimensional vector 


S Sy Mw) (4.5-21) 


Ce) 


into an m by m diagonal matrix. The operator gy l 
v 
creates an m-dimensional vector from M such that 
* 
(*) The operator yy 1 PomamCOUGUDELOUMO mGiCmArbnay agile 
lA 
operator Q (see Falkoff and Iverson (1968) pp. 3.38- 
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BM) gee cee) Lee 2.2 meee (4.5-22) 
Using (4.5-21) and (4.5-22) we can rewrite (4.5-20) as 


E[A MA‘ = B(S Q) 3 (M)) = [6, 


Ge jal 2 em 45-23) 


Proof: Letting A' = P= [Pj 5] and using the notation 


introduced in Chapter 2, we find 


m m 
(MA Nae = a Ts 1 Me = ce Ms aay (4.5-24) 
Next the product AM A' is such that 
m m m 
AMA').. = P (MA! . = : ; 
( lye nan aig! gs ee fa Bie Fak ok 


(40 5—25) 


Now we take the expected value of AM A' and consider 


thesotrr—diagonaleelements of this matrix, i.e. 


m m 
E [(A M A) i Sg oe Pee) (4.5-26) 


sky pe Bel og id 
OG equivalently 


EL(A MA‘); .] = 2 Sp in asi] =0, ij 


CARS iis) 


re 


yt * | 


ya 


pie) 


The expected value of the diagonal element of AM A'! 


is 
m 5s m m 
nPop) bee | ee aes ne + meey laa ara. m2] 
Ten eal se OM) (eat Peal ae ae RA 
k#2 
(42 5—23)) 
OG 
m 5 m m 
EF ECAR M GAs) te | Gee eee Eile. | +) ae) EF fazesmear elt. 
ii Ce. QL 12 areal Rk Lek 
k#2 
(45=29)) 
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Since Sas anda for i#% and j7#k are independent (4.5-29) 


Rk 


Cane Deureworceenmas 


2 


m 
1 — 
eee mT ee Site 


Themconclusionm(45—-25))1s derived trom. (4275-27) tand 


(A330). 


Kore chesvectorns ba — (b;) let the E[b] = b, and 


var(b) = 4 = (\;), and therefore 


Var (bebe e—weeloeb) | - = [b] E' [bl (4.5-31) 


OG 


BC (ee Ee Db je b f (4.5-32) 


=O =O 
where (A) is a diagonal matrix of the form 


PAO) We TR Ps ee (4.5-33) 


SEs pe aii, ICH Sih) 
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We now return to the problem of expressing E[x x'] 
aseestunctilon,o0L the means and variances’ of A and b. 


Usingmi4. 5-52) ands 4,.50-1/)), we obtain 


E[k ¥'] = M E[b b'] M' +M E[A M(Q(A) + bo bo) MIA‘) M! 
(4.5-34) 
or 
ELx x= M El>abilaMiee M E(A.M 30°(\) gM’ A‘) M’ 
eMeEm | AGMeb es bie 8M A] Mt” (4.5-35) 
19) =) 


We can apply the preceding lemma to U = M A(X) M' and 
V=Mb_b‘ M' and obtain 
=O =o 
Heese k eee Ome (on Se n(0)) tei S Nie ye Mae 
ry Shae 0 
(4.5-36) 


From (4.5-36) and (4.5-11) the variance-covariance matrix 


(4.5-12) becomes 


Nn 


Vacixex) —aMle [oebe ie 4 (S B10) + A(S @, 1 (v)) IM! 


= iit ey be) on | | (4.5-37) 


or 


var (¥ ¥') = M {Eb b'] - b 


+ X(S Sy 13 Mua. (4.5-38) 
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Binally eit. thes identity +(4.5-32) is used in (4.5-38) 


we obtain the result 


Var (x x') = M {M(A) + &(S 8), (0)) + eS Qo V))i Mm, 
(4.5.39) 
where 
| eT 4.5-40 
A = shy i ri ( I~ ) 
m 
O(S 8,7 (0)) = [6 , ar 055 (0) 55] ; 
(4.5-41) 
Un= iM?) aM ay 
and 
US § by: 
(4.5-42) 
WS this Te iu 
AO) =O) 


In conclusion, we now have an expression (4.5-39) for 


an 
GiemVvcianCe cOvadisgloice: Watt x On xewitchelsmstavedn in 


Geis gO see lem iow Medic wInGmVabrianices Or Avand sD. 


Also the formula for the Var (x xu) is fairly simple 


Involving On lvyeMmatiisx algebra. 


We will now consider the problem discussed by 
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estimator x POs on omic OOl en. 1S 


xX = Mb-MAMb , (4,5-43) 
where: Ag ae) — Mo and M = Mo” , and x is a linear com- 
Disniocronkon stiesnormalevari ates) insAtmerlhus Sait x = (x5), 


then Xa Ne eerie wos norma Livynduistripbuteds with 


ra m 
E Esa) =e Wid as) Mae eee lee, (4.5-44) 


(4e5—45)) 
We would obtain the same result if we used Tintner's 
estimator in place of the first order estimator. 
If as in Tintner and Raghavan (1965), we assume 
theateonlyebiscerandom jandgythat Tajj= Moe then the stochas- 


tic system may be written as 


My ne?) ' ieee, hae, es (4.5-46) 


Phosmthe tarste orden,estimatorgyx; i as 


Sie oi = he © tae wees (4.5-47) 
- o = fe) Ofic ony — 
WhiilChmns stbes Same ease thes vector: xgingt4.5-46). If be, 
where b = (b.), PoenoOmnawly Olst ri buceas tenen a Loe 
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m 
2 
IIS Sp = Ae, BIS ~5- 
[ 43 as ain [ iad (4.5-48) 
and 
ee et 
Vari xe) eee May r : (4.5-49) 
J eegyl | abs hs 


Again if we started with (4.3-5) instead of (4.5-43) then 
we could still derive the same results (4.5-48) and 
(4.5-49).. 

The only difference between Tintner's approxima- 
tion (4.3-5) and the first order approximation developed 
Rese sOCCUursawhecns Othe andsbraresrandom variables. =if 
we accept the arguments which lead to derivation of x, 
then it does not make sense to use un except when either 
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x, LOGeCOn ctl lenceymUlcCneOleneceso 1 CYS Wenlustamegeoce 


the combined influence of A and b since x, is merely a 


ie 
Pines tecOno nia atOneOmeneands) . sWhereas x includes a 
eross product term of A and b. Notice that there is no 
reason why a higher order estimator could not be developed 
Boome (Ao ea nOmces Lill abeanalySslomcOpDlied tO Obtain 


its mean vector and variance-covariance matrix. 
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Th 


Cormtiiceco ition xsof the system of stochastic linear 


equations 
(eo expe ; (aoa 2) 
with known probability density function P(A, b). It 


would be ideal to test this estimator on a "real world" 
problem. Unfortunately, because of the lack ne famili- 
arity with such stochastic models very few examples of 
them have been explored. This does not infer that such 
problems do not exist, but as was pointed out earlier, 
the solution of a stochastic model is very complicated 
and most persons would prefer to ignore the stochastic 
nature of the problem and solve the deterministic model. 
Thus to test the approximation (4.6-1) and to analyse 
the stochastic model (4.6-2), a stochastic model was 
simulated. 

There are two methods of simulating the model 
(456-7) sae licmloncemearlo method ycould be used to 
generate problems of the form (4.6-2), to solve these 
problems, and to collect the data on the system's 
Solution.) Wevcan use the data collected from these 
experiments to predict the behaviour of the model. 

The problem generation is done by selecting randomly 
with the aid of a pseudo-random number generator, a 
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density function of A and b. The pseudo-random number 
generator is described in Appendix 1, along with the 
documentation of the random number tests. The asters 
method of simulation is the method used in most of 
Tintner's work (1955, 1965) and involves approximating 
the given probability density function with a step 
function, i.e., the continuous density is replaced by 


a discrete density, 


2 
+ 

0 < Be Og = giace glad Sypocens) ie Se 
(4.6-3) 

where! the value set S = {(A,, b,)|i=l,2,...,k} is used 


to generate all possible problems. This method may be 
, ; m?2+m 
viewed as a sampling procedure in Re , where certain 
points are selected and assigned a specific weight. We 
will refer to this method of simulation as "Sampling". 
ApMOnGeECarLoOssimulation was periLormed ony the 
following example. 
Basis #1 


Solve 


(M +A) x=b , (4.6-4) 


where A is a 2 by 2 matrix of independent normal variates 


with mean zero and variance 
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FIGURE 4-1 


Average-Value Model for Basis #1 
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Vari it 0.04 ORIOL : (4.6-5) 
0.09 0.0016 


The matrix Mo Tomag em OYie2 Matic x .OLeCOnstants, 


ees ; i F (4.6-6) 


and bis a 2-dimensionalj vector, of independent normal 
variates with mean b = (15, 10) and variance Var (b)= 
(0.25, 0.36). The components of A are independent of 


thencomponentssO1e). eine average value moded, 


ihe Sy Se eee, (46-7) 


is depicted in Figure 4-1. 
PiewemusemthiemiursteordersestimnatoreiLor Basis el 


Wen ou Lnidstnes means x iS 


E(x] = [2.5 (4.6-8) 
7.5 


aN 
and the variance-covariance matrix of x is 


var (x x') eet -1.04005] . (4.6-9) 


=U L005 2.38149 


To simulate this model it was necessary to generate 
dalideiewent mcd amVvectOon, bt. irom the distributions defined 
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normal distribution which was accomplished by using the 
Box-Muller algorithm (Box and Muller (1958)). After each 
PECblemMerseqeneraten, xs 1S) TOluld, and the data are col— 
roClGwmiOrmeticencoal,;mVvarltance and covariance: Of 3." A 
two-way frequency table for x was also maintained. As 
well as solving and accumulating data for x, we also 


recorded the values of 


i . | 
= {I+ f (-1)J (wm a*)J} mM be (4.6-10) 
ie 


1X > 


a 
for i=l,2,3,4,8,16. This model was then simulated for 


1000 systems of random equations. The means, variances, 


and covariances are shown in Table 4-l. 
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Results fou Basisai ll Using 1000 Random Equations 
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Inewapproximation.0 in Table 4-l/is the ESE: x 
while the approximation * is the theoretical values 
fou ze calculated from the formulae of the previous 
section. Notice that the moments of the simulated 
first order estimator agree with the theoretical values. 
We can also see the convergence of the moments of the 
ee ordergestimator to \the momentsgoftex as’ i increases. 
The last two entries in the table will be explained in 
the next section. 

Table 4-2 contains the results of generating 20000 
SVS VensmoLerandomecquatlonsetors Basis #1. eOnly x and x 
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were generated because of the time required. 
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Results For Basis #1 Using 20000 Random Equations 
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This experiment shows even better the agreement between 
the simulated moments and the theoretical moments of sae 
Ther two-way frequency tables for x and Se using 20000 
systems of random equations are given in Tables 4-3 and 
4-4, respectively. The number in each cell represents 
Ehepnumber Ore times X or xy (eSeclemcascelay—ebe)— tai ic 


THLOetChetee Colm nes boundary spoints Of Gachycell appear 


on the appropriate axes. 
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Results For Basis #4 Using 1000 Random Equations 
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TABLE 4-6 


Results For Basis #5 Using 1000 Random Equations 


Variance 
Cov (x, ,x,) 
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Tables 4-5 and 4-6 are similar to Table 4-1 except 
thatethey icontain thercornesponding-results for the exam-— 
ples referred to as Basis #4 and Basis #5, respectively. 

Basis #4 

This example is also defined by the model (4.6-4) 


except that 


Var(A) = pees aH (4.6-11) 
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and 


M = LEB 10 : (An Ome) 2) 
0 0.0 


Basis #5 
The Basis #5 model is again the same as Basis #1 


except that 


Var(A) = O01 0.0 (4.6-13) 
bee 4 
and 
My oa ee) 0.0 “ (4.6-14) 
ee ibe 


Again these models show the trends pointed out 
earlier in Table 4-1 regarding the convergence of higher 
order estimators and the behaviour of the first order 
estimator. In summary, these results show that the first 
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and that the formulae for the moments of this estimator 
are correct. The results also indicate that a more 
accurate approximation could be given if a higher order 
estimator were fs However, as is obvious from the 
Jas easececuLron  tieacalculation Of the theoretical moments 
of higher order estimator is very difficult. Thus the 
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Tintner's Sampling technique (1955) was used for 
the Basis #1 model to compare its results with those of 
the Monte Carlo method. Discrete approximations were 
USedstOnr the Standard normal distribution. These 
apDpLOxiMNaelOnsm Involved the use of 93,055,501 /) points 
with the corresponding probabilities and these approxi- 


mations are shown in Table 4-7. 


TABLE 4-7 


Discrcteslenciuyerunct ons) Tonysampling) Technique 


SF ivewive 5 ere yahiee 
X Pr (X=x) 
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The approximations are translated to fit the appropriate 
normal distributions. Note that if there are k possible 
VoiUlesmaGQi Cac Olle m+ Me Variables (A and b) then 


2 
there are k™ we possible problems to be solved, i.e., 


83 


(2p lo oes ange WiOt9s problems for the 3, 5, and 7 

point approximations, respectively. Thus it is apparent 
that the number of possible problems is large even when 
tiemsValvlcseOremedidaksace small. Also the results of the 
sampling experiments are biased by the particular appro- 
ximation used, i.e., increasing the number of points 
sampled does not necessarily improve the results of the 
Simulation.) A comparison of the results for x, generated 
by the Monte Carlo and the Sampling techniques is shown 


in Table 4-8. 
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Comparison of Monte Carlo and Sampling 
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A) / ereeCescoatveCOnaitioneron, The First Order Estimator 


As was pointed out in Section 4.5 a necessary condi- 
tion for the existence of the first order estimator x is 


that 
Pra AM cece! ies near 1 a a. (4.7-1) 
Recall that A is a matrix of independent random variables 


ands M = Mot 1 Soma cigioor OL CONStante mel nus» for the:con— 


dition (4. /-))etospesuserul,we must first define the 


Speci f1canocmeinccmeconstderatvion. The norm in) (45/-1) 


will bepdefaneds tosbe theyHuclideansnorm, 1.e., 


m™m m ie 
Wiad) =k. i eee (4.7-2) 
i=l j=l 7 
or equivalently 
m m 
iit cee eee ce (467-3) 


i=l j=l 7J 


NotI ce the Buclvdean normesatasties the conditions (4,4=-6) 
through (4.4-10). Secondly, we need the following theorem 
which is stated here without proof (see Hogg and Craig 


(1968) p.47): 


Theorem 7 
Let u(X) be a nonnegative function of the random 
Variable. selinl wus imexasts;, then, for every positive 


CONS tallunc; 
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E[u(X) ] 


Pr (a Oe)P 24c))i< - 


(4.7-4) 


Clearly ||A ||? is a nonnegative function of the 


random variable A M, and 
m 
(A Ee = z (4.7-5) 
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Thisssubst Vel tinge nemtormuday (4. /-5)\anto the definition 


of the norm, we obtain 
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Taking the expected value of (4.7-7), we have 


2 Meee i > 52 Gas? 
OE EN ANE ee Mes OF, 7 (4.7-8) 
i=l j=l k=1 “J 
s = } eee > 
since lhe eo I = 10%, Wake Vals 9) = Oise and Beles aro! ae 16) 


ie dede eigyel 465). 

The condition that the Pr(||A M|| « 1) be near one 
is equivalent to saying that Pr(||A M|| > 1) be near zero. 
For this latter probability we can apply Theorem 7 with 


Bll elatel Wy es AEA m| |? and obtain 


thet.) 


ris 
A San aa . : 
si2 To moigoausy sv isepancort) Ss 44 iy Aj vyixrsel> é 
bas .M A sigsdigev mGGnee 
‘ —_ 
iti 
(e-F bb) 2h ee = (M A) 
a : i # ‘« {= = 
® - . Ps ; 
qatesaitgeb sid otai (7-S.6) slimmot off noijpeeedae eee 
a f 5c =3 ve ei peealy € oxld 30 7 
c m m in e 
S=f A) 1 } >= he ay 
{=; = i= f 
x0 
m o P int mn im : _ 
a rae 4 ts on - 5 
Pe 008 ake Sk i ch Lott +, - . it a 
— "ij 7 =f we & Pd ; f 
= 7 
(T-T 2) ' 
ed) aw (F=¥ ti bojosque ertd pre 7 
. : a. ni ! im c 
-_ oe *“ ; , ~ ae 
(8 ,.*) 4 of is yt 2 mY a = { 1m A} {7a 
; l=odi [=e [=r 
Al oc: : Z S 
Y = b pig 64/40 -bns ‘75> = ( fi) maV 25 by pelt 
ci.et G 7 
oe SE bis Bed 7 
ano "450.509 ({{[ > 2 Al li pbeede sxe sii | rh itz 
7 1) Ap) p34) S82 Sats KO IPAGS ene 
Tog y 
S358 toon ead (f < 


dstw’t maroon? eee fied. any MTEL Este cto stfu « 


_ 


Aaitsid6 Has 


*} [waa = 


icf 
a) 0 Dats £ 
"(nd 

7 


Fi 


_ ie 


> _ 


Un! 1s eer 


86 
2 2 , 
Dia G | AeM ie) cee>meime EE AC M| [oj ©. (4.7-9) 


This inequality can be rewritten using (4.7-8) to obtain 


aa) m m 2) 
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We can again rewrite this last relation in the form 


Brae Aa Mie eee) ee lee bey eee 
ye 


ah C4es7 ole) 


Kj ake 
Therefore, if the right-hand side of (4.7-11) is near one 


Ehnene! Gamakesmsensesto,use the first order estimator. 


The usefulness of the measure 
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can be seen by comparing the rate of convergence of the 
higher order estimators x LOwxe elnecompal tig. lables. 
with Table 4-9, we see that the value of s for Basis #1 
is 0.90600 and the convergence is very good, while the 
value of s for Basis #2 is 0.87000 and the convergence 
LseiOteassGO0dem==NOtace that the relatives frequency of 
the event ||A M// <1 -is shown at ‘the bottom of Tables 
4-1, 4-5, 4-6, and 4-9. Notice also that this relative 
frequency is always greater than the corresponding value 


Cus mas Cm liogecaucdmbyathe relation (427-11). 
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TABLE 4-9 


Results for Basis #2 Using 1000 Random Equations 


Variance 


Cou) oe) 


10.00000 
DOR Os 99 
EOe2 3056 
Oro Ose 
eS ae 
Teles SOR ys 


Lew Or 


a2 008 -26652| 1064.38) 9322.40. 


Relative Frequency of{||A M|| < 1}= 1.00000 


Pred Aen | a < > 0.87000 
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CHAPTER 5 
STOCHASTIC LINEAR PROGRAMMING 


Deal CrLOGUCtC LON 


We saw in Chapter 4 that a stochastic linear program 
is simply a classical linear program in which some elements 
are random variables. A stochastic linear programming 
model may be defined by (4.1-1), (4.1-2), (4.1-3), and 
(4,1-4).. The solution to this model will be the specifi- 
GablonsOLrethbesarstribuLion,of the obj;ectivesfunction, 
Stochastic linear programming may also be used to analyse 
the implications of the stochastic model as it is related 
to decision-making under risk (see Section 5.4). We found 
in the previous chapter that very little progress has been 
made in the analytic solution of the system of stochastic 
linear equations, and thus there has been even less success 
inethesanaivtrcesolluimonwoLethe, Stochastic linear prognam. 
Tintner (1955) suggests an analytic approach to the solu- 
C1OuROltPamsStociastlc linear prodram, .and this will be the 
LOU CEO the next section. However, as was found for a 
system of stochastic linear equations, the solution of a 
stochastic linear program is obtained by a numerical 
SepoLoximactony )l.e., an approximate Solution is found by 
simulating (or sampling) the stochastic program. 

The stochastic linear programming model has been 


applied in several areas of economic planning. Tintner 
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(1960) used this model to analyse agricultural production, 
and Tintner and Raghavan (1965) used this model ina 

| dynamic planning model for India. Charnes and Cooper 
(1964) also applied stochastic programming to critical 
path scheduling in which the task completion times were 


random. 


Deeeeniaiyt lc poOLuLLOn 


thewonly attempe.at anv analytic solution to a 
stochastic linear program was presented by Tintner (1955). 
Unfortunately, the apy success his method has had is in 
solving a single-constraint model with two decision 
variates. This lack of success is due to the method's 
complexity. Hence the discussion of this method is 
included for completeness, rather than as a suggested 
DEacCterCcalwsolUciTOn. 


Let us define the stochastic linear program as 


Zeer Tcl Some ae ge CC 22-1) 
x x 
subject to 
AS == eb Voree 2) 
and 
0 a | (5.2-3) 


Where Awdis an m by nim matrix, and b, ¢, and x are 
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dimensioned accordingly. This model is the same as the 
model introduced in Section 4.1 except that slack varia- 
Dblesiiave beensadded. | The elements A, bb, and, c¢ may be 


random, with the known joint probability density function 


Re (Mtn) m+2m+n 


Meio i, 2) aa | a eee ae (5.2-4) 


The set S will be referred to as the parameter space. 


For such a stochastic linear programming problem, 


ehenewere she Ga possible bases { BY); el cawe 
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ala) 


basis, then we can define the set 
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ponding to the k 


Sao 5S Mase thesreqicneinewnich the ee selection-.is 
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feasible, i.e., 
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Tintner claims that the regions Sir S3 


nye ClmeLia Dlcm (AyD aC) m= INiymetierencxistss ane esuch 
that 
zits = jibe z ‘K) , Kee a rcteherpan (5.2-6) 
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We can then define a region e Suchmienal 


(k) 


ue = {n | 2‘) = max 2 andes eos, (Deed) 


on 


or TY. isp the. region in which the oe selection is maximal. 


Tintner claims. that the, regions T ale eee ee Oe 110 & 
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overlap and thus the regions s.9 T. = Ur, 


Caio, chen ithe Kt 


Shee ere US 
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COBNOMsOVC Ha mee tery ce LU selection is 


k 
both feasible and optimal. 


Stochastic linear programming maps the set S into 


Re where this mapping may be a many-to-one mapping. 


K 


However, S= UL where U,_, 
k=1 : 


disjoint, and the mapping from U 


Ke uta EOLeRMuULiUa WLy 


k to Re is one-to-one. 


The transformation from U, to Re is the same problem dis- 


cussed in Chapter 4 where 


ghd coil (GN) a passe (5.2-8) 


and 
(k) 4 (Cus): (k) (5s) 


We have shown how to determine the probability density 


oe and Strom this “density, 1 wis possible 


to derive the density of Fe Thus we can derive K 


BUNCE ON Oia x 


one-to-one mappings of the mutually disjoint, exhaustive 
subsets UL Om on LO Re wealhnis again permits us to use 
the transformation of variables technique to obtain the 
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If the reader thought the transformations of the 
previous chapter were complicated, he is advised to 
ScTUuGvel 1 ninem saexamplen((1955)epp.) 205-15), of@this 
techna ques torenwmeunies i(1,02)-% Againaatishouldsbe 
stressed that while the above arguments may lead to an 
analytic solution, their complexity eliminates any value 
they might have as a practical procedure. Stochastic 
programming can best be solved by using an approximation 
obtained by simulating the model. 

However, before proceeding with a discussion on 
simulation, we will digress to present an approximation 
which has not been mentioned in the literature. It is 
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possible to-approximate the distribution z On a 


Bare CUlargbastemondgtoecalcullatenthe (meanvandivakiance 
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(Ope z (K) Why not determine the expected value of E[z ] 
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Thesprebabilityethateasparticularibasis is optimal.is 
unknown and is still an unsolved problem. However, to 
illustrate the above method, we will set i equal to 
the relative frequency of the con basis being optimal. 
We will use the figures obtained from simulating 19000 
random linear programs corresponding to Basis #1 (i.e., 
the equality of (4.6-4) is replaced by the inequality 
less-~than-or-equal=—to) > The values of Pye ke oan = 5 
are taken to be the relative frequencies of the pairs 
(x), Xo)y (Xo, X3), and (x,, Xy)y respectively, obtained 


Proms lable 5-2.) The values of sks 


eb Te Ghetel Gy (erebe, Jere, 
taken from (row 0) Tables 4-1, 4-5, and 4-6, respectively. 
‘The relative frequencies of the bases k=2, 3, and 6 are 

zero. Using these figures, we can calculate the approxi— 


Mater valuemors(>.2—-10)) and (5.2-11).. These values are 


contained in Table 5-1. 
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We could adapt this method by employing the first-order 


(k) (k) 


estimator x sigg jeWlleler=) Tepe >: elicm last rowlOL lable 


Daleshows the results of simulating 19000 linear programs. 
We can see that this method provides a crude estimate of 


the moments of z*. 


Se Om oll latLOnmOL oLOcnastiCc Programs 


Thesapproximatesdustem bution sof Zz may be found for 


Zoe eile Xam Ca x (5.3-1) 


1(5.3-2) 


= 
+ 
DP 
x 
4 
(oF 


and 


by simulating such problems according to the given proba- 
Digcyecetsatyerunction. [iis proceduremas, similar to 

the method employed in Section 4.6. We select at random 
points from a population defined by (5.2-4), generate 

the corresponding linear programs, and then solve these 
programs retaining the pertinent data. Again this simu- 
lation technique is referred to as the Monte Carlo technique 
while the Sampling method refers to sample points and their 
weights. The methodological details of these simulations 


Canewe round in Appendix 2. 
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Tomanalysesthesbenaviour of the stochastic program, 
let us consider the example, Simplex #1, in which the 
variables A, Mor b, and c of the above model are assigned 
the same values as in Basis #1 (see Section 4.6). If we 
use a Monte Carlo simulation of 19000 such problems we 
find that the average value of z is 20.05181 and the 


Vetelance OnuzZeLs: 2.0808/,. (See-Table 5-2). 
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GeneraccdsMargindgieDistribution of x3 Using» the Monte Carlo 
Method And 19000 Simplexes 
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The average-value model of Simplex #1 is shown in Figure 
5-l. The difference between the average optimal value 

of the stochastic objective function and the optimal 

value of the average-value objective function is small. 
HOWeVe TT mCoCEpDEODab minty thatvzi< 19 7orez) >) 2) 5 is appro- 
ximately equal to 0.5, i.e., the optimal value of the 
stochastic model is not very stable. In other words, 

the probability that the absolute error, |20-z|, is greater 
ehane Wels vateleastmuns. | [ney upper bound T in the frequency 
table of z, Table 5-3, refers to the upper bound of the 
frequency class, i.e., the first three frequency classes 
“eee (cope OP mel Gye) Oe) 2 ory manda 6.62516. 25) The 
second column contains the relative frequencies of the 
intervals. The third column contains the cumulative rela- 
tive frequencies of the point T, e.g., the cumulative 

Pelat tvGeerequency =o ply. 756 1S 005395 =  hescorresponaing 
information for Xp Xor X3r and xy is contained in Tables 


54 eo pO Olan e oy yarespectuively. 
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While we cannot prove that the resulting distribu- 
evoOn cole ae S normal pelt A Sei nteresting to scompare: the 
observed distribution with a normal variate whose mean 
and variance are equal to those of the observed statistic. 
This was done for Simplex #1 using 15000 random simplexes 
and the comparison is shown in Figure 5-2. It shows that 
a normal approximation for z is reasonable. 

since intnecrealways—uses) thesSampling technique 
to simulate his stochastic models, Simplex #1 was also 
simulated using the Sampling technique. The first two 
approximations given in Table 4-7 were used and the com- 
parison of the Monte Carlo and the Sampling results are 
shown in Table 5-2. The comments made earlier about the 
Sampling technique, regarding the bias and size, are 


again applicable. 


bees Prererences Functionals 


We mentioned in Chapter 4 that stochastic programming 
could be used in problems involving decision-making under 
risk, i.e., there may exist several possible probability 
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example, if we are given two possible density functions 
Simplex #1, we would want to use the model which has 
Maximum expected value of the objective function. Let 


density defined by Simplex #1 be 


Big Ave eC eG (5.4-3) 


) 


an alternative model, Simplex #2, we define a model 


which is the same as Simplex #1 except that 


Var(A) = le oe : (5.4-4) 


Or 0:5 O70,0 2 


Var (b) COMSo5, TOR INS) (Ge) 5 ES) 


It 


and denote its density by 


Ae 1g Mey) |e (5.4-6) 


Both models were simulated and the results of the two 


experiments are shown in Table 5-8. 
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TABLE 5-8 


An Alternative To Simplex #1 Generated 
From 19000 Random Simplexes 


Statistics Simplex #1 Simplex #2 


mean Z Z0F Odo 20.02208 


variance Z 2.08087 ~03641L 


0.09491 nO 2p i 
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variance xy PES E52 afOS)a bays 


variance Xo -56189 PUBS) 


cov (x), x5) 208877 -02580 


rel.freq. (x, ,x,) . -0488 sOd29 
rel.freq. (x,,x3) 5 P09) Shale) 


rel.freq. (x,,x,) 20003 s000L 
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However, if the preference functional is of the form 


MinaVan(zeleo)e | GMen ; (5.4-8) 
6 


then we would select the alternative model. 
The preference functional may even have the 
form 


iigsed he r neh, A , (5.4-9) 
where as is such that 


Pr( zee oral Ona 0.0D5 | (5.4-10) 


i.e., we choose the distribution with the greatest value 
of the lower 5% point. The relative frequencies of the 
two distributions can be compared by viewing Figures 5-3 


and 5-4. We find that we would select Simplex #2 since 
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The above illustrations have indicated how stochastic 
programming can be used in decision-making problems under 
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5.5 Active and Passive Stochastic Programming 


A variation of the stochastic linear programming 
model introduced by Tintner (1960) will be discussed 
in this section. Tintner named the model defined by 
(a a2 oh )ene (oe —2)), eGo e=3)i; eancd.e(5e2-4) tthe. passive 
stochastic linear programming model. It is passive 
in the sense that the decision variables are chosen 
after the values of the random variables are observed. 
Thus the passive model resembles the wait-and-see model 
of Chapter 2. The alternative model, referred to as 
the active stochastic linear programming model, is 


defined as 


Te Xap Coax (555-1) 


Subject to 


iceh pahe Se PGs) AU ae (5.5-2) 


and 


arn 4 (5,.9-3) 


where A, b, and c are the same as in the passive model. 
The operator is the diagonalization operator defined 
previously. The matrix U = ied Wsera Waa OVe Mela tig 13 
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ORS i5 ee, le, Sis ll ee Slt, Qe serie, ned oe) 
and 
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2 ae SL jp alae a aris (05.51) 

j=1 
The quantity Ws 2 TenRecentsoecthne proporcion, cL) tie ay 
resource, bey allocated to the qe aetivy ty, yee The 


allocation matrix, U, is assigned values in accordance 
WEEDe (5.04 manda oo )e DeLOre «the random»events occur. 
Thus the active model Asgewiilas the here-and-now model 
Ob Chap crs emia ticmactive model sthe distribution Of 
2elceoerhUnCt OlmOLreUmasewel assorted. Hence, ai a 
preference functional is to be maximized for the active 
model, it can be maximized over the set of possible 
allocatlonsmlatrices, 1Ulaethe active formulation does 
Not) permit ausimpler analytic solution. To solve an 
active model, we must again simulate the model and use 
the results to predict the "real world" phenomena. The 
paper by Tintner and Raghavan (1965) presents a Bede ly, 
good example of such a simulation. 

Some theorems on the relationship between the active 
and passive stochastic linear programs are presented by 
Senguta, Tintner, and Morrison (1963) and by Senguta, 
Tintner, Boe Millham (1963). These theorems are similar 


tomtenoccepLesentecdminsseCtion. 2.5 7and as in Section, 2.5 
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they find that the optimal value of the active model can 
never be better than the optimal value of the passive 
model (see (2.5-8)). In another paper, Tintner, Millham, 
and Senguta (1962) extend the duality theorem of deter- 
ministic linear programming to stochastic linear program- 
ming. 

Before concluding this study, we should again 
mention that although the active-passive classification 
provides anvalternacive formulation for the stochastic 
linear programming problem, it does not simplify the 
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CONCLUSION 


If a "real world" problem can be modelled as a 
linear programming model under uncertainty or a chance- 
constrained model, then we have shown (Chapters 2 and 
3) that there are acceptable techniques for solving 
these models. If on the other hand, a stochastic linear 
programming model must be used, then the easiest appro- 
Xximation to the solution is found by simulating the 
model. Although it may be theoretically possible to 
ObtraivuedleanalytlcesOLuciOn OL an, algebraic approximation 
for the decision vector in the stochastic linear program- 
ming model, it appears that the best alternative to 
analyse the problem is to simulate the stochastic model. 
However, if the problem that we wish to solve can be 
expressed as a system of stochastic linear equations 
whiCcheSatasryetoesconditions of Section 4.77, then we can 
either use the approximation proposed in Chapter 4 or 
simulate the system of stochastic linear equations. 

The first order estimator X proposed in this thesis 


is equivalent to the estimator xX, introduced by Tintner 
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(1955) except when both A and b are random variables. 
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extend, to iigher sorder approximations than is the case 
for x. Also using X as an estimator, we were able to 
Ceti Ve@iiesectiong4,/ a ysuLficlent: Condition tor x) in 
any particular stochastic model. We also verified that 
x is a good approximation for x by simulating a number 
of different systems of stochastic linear equations. 
Furthermore, this study has shown that APL in. 
its present implementation is not suitable for large 
or lengthy simulations. However, by using the APL 
random number generator in assembler form, we have 
verified that this random number generator is satis- 
PaAcCeoOuvE Om Liao LyPeCEOLusIMulalilon, 2.6.) Lhe simulations 
confirm the results of Appendix l. 
While it is hoped that this study has drawn together 
Many existing techniques and added a few of its own, 
there are still several unanswered questions. One of 
these questions is how useful are stochastic linear 
programs in "real world" situations. Although the average 
Vries OmeulesOplilalmSOLutson Of aestochaseic linear 
programming model may be very close to the OpLigial solu— 
tion of the corresponding average-value model, the optimal 
solution of the stochastic model may be very unstable. 
Hence researchers should be encouraged to incorporate 
SuOciastLc valrlattons into their models so that the merits 


Optics stocnastic lmmear programs may be assessed. 
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Probably the most interesting unanswered question 
iomuneqdetermination Ofetherprobability, densitywfunction, 


foe of the = order estimator, 


n : 7 
a ee) (MA) Mab ue (6-1) 


Once fn has been determined then it may be possible to 
show that Eid Peace ec Where fo 1Sethe sprobabi 1 ty, 
Dens wey oUnCULOnEO ma AlSO ui tf the gorobabi lity ior any 
given basis being optimal were known, then it would be 
practical to find approximations of the type mentioned 
in Chapter 5. Furthermore, we could make the simulation 
of stochastic linear programs more efficient by using 
parametric programming techniques instead of always using 
the simplex algorithm to solve each linear program 
simulated. 

As the previous pages have shown, the solution of 
the stochastic linear program is a complicated procedure, 
and until simpler procedures are found, the stochastic 
linear program will not become an important "real world" 
model. However, only the use of the stochastic linear 
program will determine its merit as a "real world" model 


and encourage further research. 
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APPENDIX 1 


TESTING OF THE APL RANDOM NUMBER GENERATOR 


Nile neroduucti1on 


The object of this appendix is to describe the 
testing of the APL random number generator in order to 
FWUSELEVeENG USCeOneARE for simulations. Thes;only other 
testing of the APL random number generator appears in 
thespaper by Lewis, Goodman, and Miller (1969) which uses 
Good"s period, Fourier-transform, and run tests on the 
numbers generated by the APL random number generator. 
Although the tests used in this appendix are different, 
the conclusion is the same that the random number generator 
is very good. The paper by Gorenstein (1967) presents 
results for the GPSS (General Purpose Simulation System) 
random number generator, while the bibliography of this 
paper and the bibliography of Hull and Dobell (1962) and 
of Jansson (1966) contain references to other random 
number generators. 

A random number generator is a device, which produces 
numbers in some nondeterministic way. The sequence of 
numbers appearing when a roulette wheel is spun a number 
of timles is random since we do not know a priori which of 


the numbers will, appear on any particular Spin. This method 
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of generation could be simulated by the computer if we 
equipped it with a mechanical roulette wheel or stored 
a table of roulette wheel spins. These two methods are 
not used because the former makes it impossible to 
repeat an experiment, while the latter requires too 
much time and storage. 

Instead, a recursive mathematical method is used 
to generate a sequence of numbers with certain properties. 
Since any number in the sequence may be determined pro- 
vided its predecessor and the method are known, the numbers 
abemcalledsspseude-rendom numbers = lhe particular’ recur 
sive method used in APL is due to Lehmer (see Hutchinson 
(1966)) and may be defined as follows. 


For parameters k and x, we may define a sequence of 


1 


numbers x oe ~, as 


ee! 


= 18 xe (mod p) A ig atl a Array 5 Oils Sah) 


OGeLneAPLenotation 
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Sucnpaegenerator isereferred) tocas, aimultiplicativesy cons 
gruential pseudo-random number generator. The parameter 
Kee SREhesmUbtiplaer, xy is the seed, and p is the modulus. 


Thescharacterirstics=of the resulting stringfof numbers,is 
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dependent upon the particular set of parameters chosen 
as is shown by the following example from Hillier and 
Lieberman (1967, p.446). 


igs: X= See) Gaps) LO, and 


x 16 qi x (mod p) gt a Re OF 


then the resulting string of numbers is 


Notice this particular combination of parameters yields 

a sequence which repeats itself every four digits. It 
would seem reasonable, that among other properties of the 
generator, we would want the cycle length as long as 
possible. The cycle length is determined by the values 
Of kK and peaccondingeto ssome number theoretic theorems. 
However, these theorems do not yield a unique set of 
values for k cee ay but rather they suggest optimum pro- 
DebliecssObekeand pw ePirst peshould betprime land fas barge 
as possible. Since the word length of the APL interpreter 
for the IBM System/ 360 is thirty-two bits, the maximum 
Value storap 1s Das - 1, which conveniently happens to be 
Disttich we yliecimemeshOuld $be;, ajpositive primitive rootio ip 
and should be close eye) Lewis, Goodman, and Miller 


(1969) show that ee POCOverS £SUGI IomVy cae. for eke erhis 
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the APL random number generator is 


5 Syl 
ee ee di Lenters = Aley al Ay eG op (Arye hee) 


OF 


JG 2 MET Pe DS HIG hae ee ae 
This particular generator has cycle length of ape = Yh. 
The seed in a clear workspace (the workspace obtained 
bys the Wiel FAR command ) eee same as the multiplier, 
16807. APL then saves the last random number generated 
and uses this number as aK (OX eb Or the nex srandom 
number generation. A )SAVE command will store the last 
random number generated for the seed when the workspace 
ismeceloadeds @)0AD. 7). 

We now have obtained a random generator (Al.1-2) 
Witchec Maximumecvyole length. 7What other properties do 
we want of a generator? First, we want a simple, short, 
and fast generating procedure, which obviously (Al.1-2) 
is. The paper by Lewis, Goodman, and Miller (1969) 
Meneions an Uppexr bound of approximately 31.2 microseconds 
per random number call in FORTRAN for a System / 360, 
Model 67. Secondly, the particular method of generation 
must be reproducible, a necessary condition if we are to 


be able to repeat simulation experiments. We can do this 


1S 


with the method (Al.1-2). Lastly, the numbers generated 
must be statistically acceptable,’ i.e., they must: appear 

to be random numbers. The particular statistical proper- 
ties which are most important will depend upon the par- 
ticular application employing the random numbers. Usually 
the random numbers are used to simulate a given distribu- 
tion. The question then is how well do these random numbers 


Simulate the distribution in question. 


Rie 2 ee new APL RandomaNumber Generator 


Before considering the specific tests, let us first 
look at the testing procedure used in APL. The code used 
in’ the interpreter to implement the generator just des- 
Gribed is not accessible to the APL user. ©The user has 
access only to the operator ? which he may ‘use either as 
a monadic or clyadic operator. We are interested only in 
Eheemoneavem: @. em lLeeke1Sea=pOSsSltive=integer scalar, then 
?K produces one number drawn at random from the set 
Lo ee ea el ne) eee tThussthesanumbers generated 
Dyeerk Shioulanbesvallcsurrom=the diserete™ uniiorm dis= 
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It is the acceptability of the random number generator 
in this form in which we are interested, since any simula- 
tion in APL must ultimately be defined in terms of the 
MOnad1cey se Thess value of K used was 10, and the number 
2K was added to ‘1 in order to produce a value in the 
So Ome en nie et neepart culo esrorn-oieet ne 
generator was 1+ ? 10. The starting value chosen as 
the seed was 16807. In other words, the tests were performed 
on the random numbers that would be found by working in a 
clear workspace. The APL function written to test the random 
numbers was TEST with the header W TEST M where N and 
M are positive integer Sscalamusamelhe Lirsting@otethi Ss func— 
tionsmay be-found at the end of this appendix. This function 
is self-contained except for the gobal variable PRINT (a 
positive integer scalar) and the subroutine RWD (a rounding 
hOuUgIie pESee eodams a (L967 )iepDso—/).. Une LUunCtIOnstesesa lea 
random numbers as well as each subset of PRINT x M random 
numbers where M is less than or equal to 1000 (because of 
Pinmecemscace mandua multiple of 4° (because of the Poker 
Desc) Nemrunch On tceypes a repore Oneeach Subset (of 
PRINT x M random numbers as well as a summary on all WN x M 
random numbers. The particular values of (N, M, PRINT) 

Meo Cmve Gem Ome UUPme ln) em 0), LOO0); >) meeatnd ( 100," 1000;, 310). 
Onlvethemlast results of each run willvbe reported here, 
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Timing Estimates For The Function 
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Ree LeSCeIDelOnmOreotatistical Tests 


In this section we shall discuss the tests that 
were used to test the monadic operator ? . The six 
tests were the frequency test, odd and even test, poker 
test, dap test mcoupon stest, and serials test. “Each will 


be described briefly 


a) Frequency Test 

The assumption that the numbers produced by 1 + ? 10 
are uniformly distributed can be tested statistically by 
uSing a frequency test. The number of occurrences of 
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where = Pope NemnUnlet Of a .Seinethes vector Die This 
result is contained in the first frequency table typed 
Dyn Oo MeaWwitcnmcOnatatcnor dd 10sby 4etablerinswhicem the 
PEnotecolunn comeainomene digit, thessecond column, the 
frequency F, the third column the relative frequency 
(£./M, alee gl (Orme Mj, wan stiesLouUrtoe column tne 
cumulative relative frequency. The results of this test 
for the three runs are shown in Tables Al-2, Al-3, and 
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TABLE VAL =3 


PrequencyelablesProducedsBy eeliol.. 


For 50000 Random Numbers 


Relative Cumulative Relative 
Prequency, 
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Theoretically, the numbers generated should have 


the probability density function 


Bi (X=30) 0) ee, Kooi A eeev.c es ee ne, (Al 3-3) 


= ; otherwise. 


If the numbers are random then they should have the above 
DpEecbabriitvy Censity function and thus it 1s possible to 
check the moments of this distribution against the appro- 


priate sample statistics. The first moment or mean is 
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Clase CONS stcrom onlveone point. A statistical test of 
the sample mean can be developed since the sampling dis- 
YOU CT On EO ts meus approximately normal because the sample 
size is large and the population variance is finite. 

Thus m can be considered as a value of a normal variate 
with mean uw = 4.5 and standard deviation o//Nx M = 


2.87228//N XM . One such test for the sample mean could 
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TABLE Al-5 


Conftyvdence intervals Kor The Sample 
Mean In The PRequencyetest 


Confidence Number of Confidence Sample 


Level Random Numbers Intervals Mean 
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Wy = Elx”] = 2 x7p(x) = 202.5 (7s 329) 


or 
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Piecm Values sosn thesabove «statistics for each of the three 


TUNSTOL SLES are recoraed in Table Al-—6. 


TABLE Al-6 


Brequency lestestatistics 


Observed Value 
Theoretical 
Scaersercs x] aula 100000 
Random 
Nos. 
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Recall that Tables Al-2, Al-3, and Al-4 contain the 
observed jirequency of Gach digit, and notice that by using 
(A1.3-3) and the number of points generated we can calculate 
thesexpected frequency of eachedigit. Using the observed 
and expected frequencies, we can apply the chi-square test 
to determine if the numbers observed come from the population 
defined By etl oO ReENCECCEanewa StatisuLc ve CCH i Ssuch 
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following) TableweAl-7/ shows the results ‘obtained by TEST. 


Notice the only questionable result is the last run. 


TABLE Al-7 


Chisoquare StaclouLle Forrrne Beequencvelese 


Number of Digits Critical Value} Chi-Square Statistic 
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b) Odd and Even Test 

If the numbers generated are random, then there is 
an equal chance of generating an even or an odd digit. 
In this test the number of odd digits 9, and the number 
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We can then postulate that if the numbers in D are randon, 
EhonwOnetnesaverage mM + 2. digits should be odd and M = 2 
digits should be even. Thus we can again test this hypo- 
thesis using the chi-square test where the statistic es 


(Comme secefined eas 


e712 (e 'u/2)? 4 (@ = /2)7} Gn 17) 


or 


Corrifame ane ND meee a Sa )8 Ds), teeGOH byMint Die tx. 2%. 


This Statistic also has an asymptotical chi-square distri-— 
bution with one degree of freedom. Thus one test of the 
above hypothesis using a confidence level of 90% is to 
accept the numbers as random if va < 2.71, otherwise reject 
them as not being random. The uc results of this test are 


shown in Table A1l-8 and are all acceptable at the 90% level. 
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Une output Ofethis test consists of the observed 
relative frequency of odd and even digits along with the 
value o£ the appropriate chi-square statistic. This 
output is generated by statements 67 through 72 of 
Diol. 

c) Poker Test 

Tia Ssmeect me) seecOnCermnecde wet thesdistripution) OL 
certain quadruples of these generated numbers, and hence 
is given the name "poker test". We restructure D into 
M + 4 sets of quadruples and then test these quadruples. 

We are interested in particular combinations that may 
occur in these quadruples. Specifically, we want to know 
in how many of the quadruples there are all different 

Gig eespmaOneMDaliwOmesitibear digits, two Dalvs Of Simular 
digits, three digits alike, and all four digits alike. 
BOmeoimolvcrly ,eLctuuseCcally the event of all ditterent 
digits event one, the event of one pair of digits event 
EwO, mands SO. On. 

The major problem in this test is finding an algorithm 
which will map the quadruples into the specific combinations 
using as few comparisons or tests Pence! This was 
necessary in order to economize on time and space. The 
algorithm used makes six comparisons which, if we represent 
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Thus we construct a vector vey Sukelat qelekehe 
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Hence, we have transformed the M + 4 quadruples into a 


Moz 4 component vector U such that 


6 
es (u; ) where Ure 2 Waa ate etna hoa ea! (Al .3-21) 
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which is approximately a chi-square variate with 4 degrees 
Ofeirecaom., ThUS we Can uSe this statistic to test the 
hypothesis, that the numbers are random, by accepting the 
hypothesis if 2 < 7.78 (at the 90% level) and rejecting 
it otherwise. The results produced were all acceptable 


and are shown in Table Al1-9, 


TABLE A1-9 


Chi-Square Statistic For The Poker Test 


Number of Digits Critical Value Chie Souar cma al oC 


generated . 


10000 
50000 


100000 


THeTOUGPULCNOL thls eLeSsSt COnSISts of a table of the 
SeeCliverCVeNuS eLicir tmeOoLletical probabilities, the 


observed frequencies, and the observed relative frequencies 
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and@also.the valuc lot. the ra Chiesiuaserstatist Cm he 
accumulation of these data is done by TEST in statements 
10 through 13 while the report is generated by statements 
ies) Wioawaiutete) (335) 

d) Gap Test 

Not only is the number of occurrences of any digit 
important, but in some simulations the number of digits 
between successive occurrences of this digit is also 
important. The number of digits between successive 
occurrencessoOf a pareicular digit plus one) is defined 
Combes tse dapeengcl wee rOnr example, 16 the string of random 


numbers is 


Smo cece BOM O24) > O60 OFS 
8 7)» waa 3) 


Ene emc amen GensmrOnethie digit os) dbe.0 m2, )y. ands >. 
Sincemelem=pLobabtintysoL,oCCuLrence On any. parulcular 
digGvlemcm0.J7 sacsuming the numbers are random, the 


DeObability Oreawap length xvis 


“(sb heel ea bree 2S er Wea oS) 


The distribution of the gap lengths is a geometric distri- 
DUET OnEWicdmaepDrobability Or SucCcCeSS or pe 0.1 such that 


theameaneis il=10,and the variance is oe 008 
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Given a string D = (d.), ves eM Ole rT ancom 
numbers, how do we find the gap lengths for each digit 


iyoi=Vplpe ee roe MPirste, construct a@vector es Sucnietnat 
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each j such that 
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we are able to obtain the observed frequencies and statis- 
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The function TEST produces two relative frequency 
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contains the theoretical probability, and the remaining 
five columns contain the observed relative frequencies 
Otethesdigutsce lato 4einathe first, table andthe «digits 
DecCOMJEUn tNemsecongatable, The function also gives the 
number of gaps, mean, variance, standard deviation, 
maximum gap length, probability of the maximum gap 
lengthwandether chi-square Statistic for each digit. 

The chi-square statistic is calculated on the 
basis of the lengths 1 to 23-ignoring any gap length 
GVeavepachiane.o mee becouse whesamount Of information ‘thar 
can be stored is limited, the frequency table is truncated 
gee 2oesetnesgap lengths greater, than 23 could have been 
pooled, but this would only permit another degree of 
freedom and at the same time introduce pooling errors. 
Hence there is not much to choose between the pooling 
method and the one employed here. Since the chi-square 
Stabistevcetsebasedson 25 Classes, ltyhas' 22 degrees of 
freedom. If the hypothesis, that the observed data 
approximate the geometric distribution between 1 and 
Dowel omeOulCracCepted che oe value should be less than 
30.8, otherwise the hypothesis is to be rejected. The 
values of the chi-square statistics for each digit are 
shown in Table Al-10. Notice there are two rejections 
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TABLE Al-10 


Chaaoduare Closwesk heakem iene: Ube Gap Test 


Crvoioal 


henna Value 


O 


The gap lengths are 


Sg ts) lajahys 
SOR S ar 
5 00 14. 
30.8 205 
3050 JETS 
S03 Heke 
307.3 dbp 
5 Uno Zia 
30Rs Ia 
30%8 PAN 


Number of Digits Generated 


10000 50000 100000 


vl 2072 24.0 
4 IMEI MS 53.0 
6 J@)o Ah Ward 
5 Zone SHES TS 
3 1324 ions 7 
5 PASS Se) 20.9 
A 17.4 LORS 
i Zon Askew 
i 28.4 YASS sah 
3 220 2i2 


calculated by statements 14 


tome Of ThST who lesthesreport is*generated by statements 


Bs; vets) ALO 


e) Coupon Collector 


Test 


Ties COULCGHECOMleCCLOM Lest iCalCculates  theslength 


of the string of random numbers which is necessary in 


CLIecte COCCHI Cen O1gluSeto AppeaGein Lhe string. Random 


numbersearercol lected untia, all=ten different digits have 
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been selected. The theoretical distribution of these 


set sizes, x, is 


8 
ie) = a (-1)* (oe) aean GIN 


andgnaseanmeansote29.2908and varrance of 9125.687. 

Here again the main problem is to find a satisfactory 
algorithm which wulletind the coupon set size in a<string 
D of random numbers. The method used is based on the recur- 


Sive formula 
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the vector Wea The vector as we), Wee bo! 
is defined such that we PS the index of the: first ocecur— 


TeEneer Ong inthe vector 


= (d_ pt he eaop Chale s (AS 30) 


sy can no. longer be constructed when 
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[Helmet ucmVvee toGgecontains ther coupon Set sizes. % as used 
to obtain the frequency table and the statistics. 

DiCmLestheosOperivEss testers a 2omby 4etablerine which 
the Pacst column contains set sizes, the second column 
contains the observed frequencies, the third column contains 
the observed relative frequencies, and the fourth column 
contains the theoretical probabilities. Also shown are the 
number of completed sets, the mean set size, variance, 
standard deviation, maximum set size, probability of the 
Maximum set size, and the value of the chi-square statistic. 
be is *necessary ins TEST to set the Pr’ (X=x) = O°-for x > 80. 
The above results are generated by statements 110 to 122 of 
TUS 

Thescim-squahkems catistiC 1SeCa CUulated from etne 
PrUniceatedearmser DOU tLOnes(xo=—6L0, lle... , 834) .. However, tie 
use of a chi-square test here must be viewed very critically. 


Eiem ced sOnmerOre ti tom homt leat ie where 


ee } (observed tfrequency—expected frequency) ” 


expected frequency 


X 


(Adio 5 2) 


Only tends to. a Chi-square distribution. The rate at which 
x? tends to a chi-square variate depends on the theoretical 
distribution involved as well as on the number of points 


SosecuyCOmebiiceOllgigial Gcictrioutcions (AIS3-28) mustebe well 
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behaved. Many rules regarding the theoretical distribution 
have been evolved by Cochran((1952)and (1954)) and Keeping 
GE2G2)= ehOm thersdistcibution (AL. 3-28), a unimodal dis-— 
tribution with well behaved tails, Cochran suggests that 


each frequency class has. an expected frequency of at least 
4 


v 


One -ws BlnetLhescouponed istic. bution ePr «(X10 )ie= 93.6, %.10 > 
and the expected number of completed sets is M + 29.290. 

If M is 10000 or 50000 the expected number of completed 
S@useuswo4 sore /O/perespectively. Both of these values 
yield an expected frequency of less than one for the 
frequency class X=10. Almost 100000 random numbers are 
necessary to satisfy the above criterion so some compromise 
between the theoretical and practical value of wile test 
MUStebe mage. lhe critical value of the chi-square Statis-— 
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Number of Digits Generated : N x M 


Theore- 
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£f) Serial Test 

If the numbers generated are random then the 
successor of a particular number is as likely to be one 
diguiemassanOtier, .1,e.8,8 the Successor Of a particular 
digit may be any one of the ten digits - all with equal 
MeOod Die CV ein mOrC mE cCOmMakCuCitseteSt, ati rst construct 
a sequence of pairs from the string of random numbers D 


by forming the pairs 
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ence then calculate the following chi-square statistics 


aah (Ale 30) 
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CHT At ee (t/t) ee) Ms 10) e 2). xX 108s 
and 

i= 2 ‘ ae = aaa (Al. 3-37) 
Or 


CHEM. (Ce GC r= @ Mes BOO i x 100 24 My 


where RF will have 9 degrees of freedom while x5 has 99 
degrees of freedom. It can be. shown that x3 = x - x5 

is asymptotically chi-square with (10) - 10 = 90 degrees 
of freedom and x4 = x5 - 2x5 is also asymptotically ehi-— 
square with @ne ye = 81 degrees of freedom (see Hull and 
Dobelaig €b962)))i- 

The calculations for this test are done in statements 
33eto ee and sthe report which consists of the two-way fre- 
quencystable Heand the form values of tthe :chil-square 
Sta tct bOceeuSmOLCDaiGeO ein statementsral 25: toj136) of vTESL. 
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TABLE Al-12 


Chi-Square Statistics For The Serial Test 


Number of Digits Generated 
Critical Value 


Chi-square Co 
Stati1stics (90% level) 


10000 50000 100000 


Shey CLO ES 29 
10 27 187.6 TOO. G16 
9354066 85.646 


85.4416 I UCIO TES 


Nissan Conclusions 


The following pages contain the listing of TEST and 
Gime xa lLoeO ee LiemrinctronussOoutput fLormat., Lhe, output Lormar 
was "cooked" slightly to produce a suitable output (eta tO 
Biteanieos bye lepage)™. 

Generally, the numbers in the vector Sipe 10) (for 
M<« 10000 or 50000 or 100000) appear to have random number 
properties. However, it should be noted that the "randomness" 
Otsanyepalcicular set Of such numbers cannot be guaranteed. 
One set of numbers may excel in the frequency test but may 
farle thesgap test, and thus the purpose of the simulation 


will determine if the numbers are usable. Some numbers are 
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very poor random numbers and should not be used in any 
simulation experiment. The only way to be positive, 
that the random numbers are in fact random, is to test 
them, and therefore all simulation experiments should 


contain an analysis of the random numbers used. 
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APPENDIX 2 
THE MONTE CARLO SIMULATION PROGRAM 


In this Appendix we shall briefly discuss the 
Monte Carlo simulation of stochastic linear programs 
and describe how the following program may be used to 
perform this simulation. The Monte Carlo simulation 
program collects data on the solution of randomly 
generated linear programs. Using the data obtained 
from the simulation, we hope to be able to predict the 
behaviour of the stochastic linear programming model 
(see Sections 4.6 and 5.3). 
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program can be expressed in the above form, where Mo oF 
bor v, and c are known constants, then the following 
program can be used to perform the Monte Carlo simulation. 

The following table, Table A2-1, contains the input 
or data set specifications for+the Monte Carlo simulation 
program. 


TABLE A2-1 


Detapocteopec i ucatVonsHovs Ine sMonte Camo Routine 


Variables| Columns] Type Comment 


Integer|no. of rows in ie 
Integer|no.of columns in Mo 
[NRxNC elements of M_in 


+8 | row order and°Fl0.5 
format 


[NRXNC elements of S in 
Tes | row Order and F10.5 


FOrmMat 


elements of b_ in 
EALOm eG LNans 


elements of v in 
El0 25> format 


elements of 
EL0253f0rmat 


NUMBER theeno. OL) linear 
programs™ to" be 
generated 


THe sopeLacOtrmieis the Greatest integer Operator. 
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This data set must define a stochastic linear program 
of the form outlined above. 

After the program has read in the data set, it 
initializes the frequency tables for z* and Xie 
i=l,2,..., mtn. These frequency tables depend upon 
the particular model being simulated, and therefore 
must be defined in relation to the data set. Using 
the subroutines BOX and RANDOM, the program generates 
random points and creates the corresponding linear 
programs. Statistics are kept for random points which 
are reported at the end of the simulation, so that the 
randomness of the points may be checked. The simplex 
algorithm (SIMPLX) is used to solve the linear programs 
or to discover cycling or the existence of unbounded 
solutions. If there exists an optimal feasible solution 
COMCNCm@ LI Neat pLOgbal, themetne Stdbistics ton 2* sand |x 
are collected. The simulation continues until the 
required number (NUMBER) of linear programs have been 
generated. 

Once the required number of linear programs has 
been generated, the output routine begins. The output 
from this Monte Carlo simulation begins with the average- 
value linear program corresponding to (A2-1), (A2-2), 
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subject to 
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and 
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(AZ=6) 


Next the program reports the number of linear programs 
generated, the mean and variance of z*, the mean of 


x. 
eEU 


Lote. pet eendetne, vabiance-covariance matrix 
of x. A frequency table is printed containing the 
frequency of the pair of optimal basic variables 

(x5, ae: i7#j}. Frequency tables for z* and Xey 
i=l1,2,..., mtn, are also given and are of the same 

WOR BES WEIS, Sos, Bee A Bessy, Glee! Sih 4s ihe oeleKe 
program reports the mean and variance-covariance matrix 
of the generated random points. 

The following Monte Carlo simulation program 
requires the FORTRAN subroutines: TABLE, SIMPLX, and 
BOX. The subroutine TABLE is a frequency tabling routine 
which uses a binary search technique. The subroutine 
SIMPLX is the FORTRAN implementation of the simplex 
algorithm. —1heybox-Muller algorithm 1s used in BOX to 
transform uniform random variates into standard normal 
variates. This subroutine requires the uniform random 
number generator, RANDOM, which is the assembly language 


version of the APL random number generator. 
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This Monte Carlo simulation program can be used 
as a basis for Tintner's Sampling method where only 
the method of generating and weighting the points needs 
eOnbe: Changed. The Monte Carlo program can also be 
modified to either simulate or sample systems of linear 
equations. The following pages contain the main line 
i ong OretloeMonre. Camo wprogram, and the el1stingsof 
the subroutines RANDOM, SIMPLX, TABLE, and BOX. The 
results obtained by using the Monte Carlo program to 


simulate Simplex #2 follow the program listings. 
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SIMULATIGN GPR MONTE CARLO 


TECHNIQUE FOR SOLVING 
STOCHASTIC LINEAR PROGRAMS 
G. LINKS 


(enue) (ioe) (onic 


REAL*8 TAIL(5364,2)/640*0.0D0/y MEAN/O.ODO/ sMX(4)/4*020D0/1 
IM2X1494) 2X04) 2D7T(395),EA(252),VA(2,2),EB(2),VB(2),C(2),D0C(5), 
2DX(2)3R06)_PTS(6)/6*0.0D0/, VP TS (6,6) /36*0.0D0/, SUM( 64) yM2 

INTEGER FREQ(4,4)/16*0/,INDEX(2),INT/16807/, KOUNT/0/,KORE/O/ 

DATA M2X/16*0.0D0/,M2/0.0D0/ 


READ DATA FOR PROBLEM 


Aan 


READ(5,200) NRyNC 
READ(5,201) ((EAC(T, J), J=1 »NC) »T=1,NR) 
READ(5S, 201) ((VAC1T 25) 2J=19NC),I= ENE 
READ(5,201)(EB(1),T=1,NR) 
READ(5,201)(VB(1),T=1,NR) 
READ(5,201)(C(1),1=1,NC) 
READ(5,202) NUMBER 

200 FORMAT(214) 

201 FORMAT(8F10.5) 

202 FORMAT(110) 


f 
Cees ECR §S FOP 
C 
DO 3 1=1,54 
TAIL €(13141)=16.0D0+DFLOAT(I-1)*0.12500 
A Tietea its 1 = 0-2 O00) SDEU OAT C1=19 0.052500 
TAIL (€3,151)=8-0D0+DFLOAT(I-1)*0.0625D0 
TAIL (43131)=0.5D0+DFLOAT (1-1) *0.125D0 
3 TALE G5 ly 1) =(—2 -OD0) DEL DAT (1—-1)40. 062500 
5 KOUNT=KOUNT+1 
IF(KOQUNT.GT.~NUMBER? GO TO 100 
C 
C GENERATION OF POINT 
& 


CALL BOX( (NR#NR*NC),R,INT) 
DO 75 {[=1,NR 
DO 75 J=1,NC 
K=J+ C1-1)*NR 
: R(KJ=EACT,J)#R(K)*DSORT(VA(T,J)) 
LS NT( 1, J) =RUK) 
DG 76 T=1,NR 
K=NR*NC+I1 
R(KJ=EB( LI) 4+R(K)*DSORT(VB(T)) 
76 DI(IsNR#NC+1l J=R(K) 
K=NR+NR*XNC 
DB 77 I=1,K 
PPS (IN=PIS(I)+R(1) 
DO 77 J=lyK 
crys VPTS(T,JI=VPTSCI, JI +R( 1) *RCJ) 
G 
C SET UP NEW PROBLEM 


ond 
= 
= 
t 

p< 


\0 10 of ee\t iA™ \ Pe ff) VI en 19 CHrODALTS 2 


K=NR+NC+1 169 
DO 78 I=1,K 
DC(1)=0.0D0 
PEt leiGESNG) OCI J=CAal) 
78 CONTINUE 
DO 79 I=1,NR 
79 INDEX{(I)=NC+I 
DO 80 I=l,NR 
DO 80 J=1,NR 
DT(I»,NC+J)=0.000 
IF(J.EQ.1) DT{1I,NC+J)=1.0D0 
80 CONTINUE 
MR=NR41 
MC=NR+NC#1 


CaepoULVE; NEW) PROBLEM 


CALL SIMPLX{MR,MC,DT,DC, INDEX,L,ZMAX,DX) 
PEtCESEO. 1) GO 10450 
KORE=KORE+1 
LECES ED. Cie GU Oe25 
WRITE(6,900) KOUNT 
JOOP EORMAT Uta. §PREOBEEMANO. °.110,° 1S CYCLING?) 
GOatG=s5 
als) WRITE{6,901) KOUNT 
901 FORMAT(* *,'PROBLEM NG. *,110.' HAS AN UNBGQUNDED SGLUTION') 


bHad Oy 5 

C 
CenbATASCOLLECS TON 
C 


50 MEAN=MEAN+ZMAX 
M2=M2+ZMAX*ZMAX 
K=NR+tFNC 
DO 51 [=1,K 
51 X({1I)=0.0D0 
DO 52 Y=1,NR 
ae X(INDEX( 1) )=DX(T) 
FREQLINDEX(1)-,1INDEX(2 ) J=EFREQTINDEX(1), INDEX(2))+1 
DG 53 f=1,K 
MXC T)=MX(1)4+X(1) 
53 M2X{1TsJ)=M2X(1_9 JS) 4+X(1)*X(5) 
C 
CepiABi Ea DATA 
C 
CALL TABLE(TAIL,64,ZMAX,1-0D0,1) 
DGY6Or1=2,5 
CALL, TABLE(TAIL, 64yX{1-1),1.0D0, 1) 
60 CONT TINUE 


GOetoes 
C 
C OUTPUT ROUTINE 
C 


190 K=NR+NC 


cs 
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: YM= thd 
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30 


902 


EARS 


904 


905 


906 


907 
908 
909 
eek 
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NUMBE R=NUMBER-KORE 170 
ME AN=MEAN/DFLOAT (NUMBER) 
M2=(M2/DFLOAT (NUMBER) )—-MEAN*MEAN 

DGelOl 1=1,k 

MX( 1) =MX(1) /DFLOAT( NUMBER) 

DO 102 J=1,K 

DO 102 J=1,K 

M2X( 1, J)=(M2X013 J) /DFLOAT (NUMBER) )—-MX( 1) *MX{J) 

DO 30 I=1,3 

K=I +1 

DO 30 J=Ky4 

FREQ(1,J-1L)=FREO(1,J) +FREQ(J,1) 

WRITE(6;902) 

BORMAT ELI. «MONTE CARLO METHOD FOR, SOLVING!/ 

140%, {THE SSTOCHASTIC LINEAR PROGRAM :!) 

WRIT (6,963) (C(1I)s1,f=1,NC) 

Ber AN GG so ks OMAXIMIZE Ze) WHERE. *7/90',5X%5°Z2Z = *,F12.5;" X*eIl, 
Ares OPPS? SO ey 

WRITE (6,904) (CEACT oJ) ody J=1LyNC), EBL I) y1=1,NR) 
COPUAMULON, 5 XUSUBAEIAls TOU/ (L0855X 14k l205~ 5 X's 11," + *)s 
PEs et eb hs! pols 9612 65):} 

WRITE(6,905) NUMBER 

FORMAT{*0',*AND USING ',110,* POINTS") 

WRITE(6,906) MEAN ,M2 

mien bi OS MAPPER OX. VALUE OF E(7) IS) *,F12.5/'0',* APPROX. VALUE 
PRG VAR CZ) YS4>5 F112. 5) 

WRITE(6,907)01,MX(1),1=1,4) 

Soe AnteOest Are ROMs VALUE OF SE(X4. 40 lytje1S*aF leo) 

WRITE (6,908) ((M2X01,J5),J=154) »1=1,4) 

Bee MAM OM Nt ROk seen COVARIANCE MATRIX 1S */ 1205 40E1 225, 24)) 
WRITE(6,909) 

BORA Os SEREQUENCY OF (PAIRS MATRIX!) 
WRITE(6,910)(J5 J=294) 

i Aes Oo et hie 1 Xie LL) ) 

WRITE (6,911) 01, (FREQ(1,J) ,J=1,3)3,1=1,3) 

PURMAT UO Salil s30S%s 1L0)/'O"',11, LOX, 21 5X TPO) / PO LS 5oXs 110) 
HORS 5 = 

TAIL (591.2) =TAIL(J,1,2)/DFLOAT (NUMBER) 

SUM(1)=TAIL( J31s2) 

DO 31 I1=2,64 

TAIL (Js192)=TAIL(UJS,1,2) /OFLOAT(NUMBER) 
SUMUTIJ=SUM{I—DieTAGL( 3,152) 

Het Qanen i) “Gia, Wah ew 

WRITE(6;912) 

FORMAT(*1L*',*THE APPROX. DISTRIBUTION OF Z*) 

GO 10 33 

K=J-1 

WRITE(6,913) K 

FORMAT('1",*'THE APPROX. MARGINAL DISTRIBUTION OF X?#,11) 
CONTINUE . 

WRITE(6,914) 

EY Re hal ( 0s SD Xe! lop mel REO ey (Xs TCRREIEREO S's lL OXs U1 5 LOX, 
Pre EEO ee Ma Ge Re FREQ 2! } 

NO 24 I[=1,32 

K=32+I 
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35 
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WRITE (6s91L5) 0 TAIL( Jo TP ,b) sL=ly2),SUM(1T) s{ TAILLE Ss KyN) yN=1 22), SUM{ 
BORMAUC ME MsSLEI 2s 5aeX)s LOX, 3hEL2+552X)) 

CONTINUE 

WRITE(6,916) 

BorAT (SOF / 7/108 87 2 tHeE UPPER BOUND GE VAN OPEN=-CLOSED 


TARUNTERVAL?) 


WRITE(6,917) 


GOGRMAI(*0%> "ROSFREQ. ser EAlIVESEPREQUENCYSORS THE SINTERVAL (12) 
WRITE(6,918) 


BOD Ae Guge Gens care. > CUMULATIVE RELATIVE FREQUENCY UP TO & $ 


TPeENeHRUDING THE INTERVAL |!) 


CONTINUE 
KOUNT=KOUNT-1 
K=NR+NR*NC 
DO 365 I=1,K 
PTS(I)=PTS(1T)/DFLOAT(KQUNT ) 
DO 37 {1=1,K 
DO 37 J=1,K 
Ved Cle =(VPTS{ 15) /DELOAT(KDUNT 9)—-PIS (1) *PTS{ J) 
IRITE(6;919) 
FORMAT(!1','SUMMARY GF THE RANDOM NUMBERS'///) 
WRITE (6,920) (T,PTS{1I),1=1,) 
GORA MIMOSA ae HEMMEAN VALUESOF- Rltsi2,t30 1S '5E10.3) 
WRITE(6,921)I 0(VPTS(1,J),J=1-K),7T=1,K) 
PORMAT(4—*%,* THE VARIANCE-COVARTANCE. MATRIX 15'/{?0',6(E10.3,2X) 
STOP 
END 
SOURCE STATEMENT 
RANDOM Colt 
USING *,15 INITIAL LINKAGE 
STM Zap cols 
LM Zoo) UGADSAODRESS Eo ee VARA IE 
u 59 A COMPUTE NEXT INTEGER 
N 4,0(2) Als DISH AXITI A MGD Pp) 
D 4,P 
ST 4,0(2) 
SRL 4,1 COMPUTE NEXT REAL 
A 4, CHAR 
Si 4,0(3) 
LM pL be He Oe TERMINAL LINKAGE 
MVI bP 7 OU RETURN 
BR 14 
CHAR Ot Pr OW Bie le 7 at GONS TANTS soues ARs GL oe 
ch RE melos0 7 SO A LST UN DOUBLE WORD 
DC F'214 7483647! BOUNDARY. MAKES LM 


END ’ INSTRUCTION FASTER. 
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SUBROUTINE SIMPLX(M1L,N1,As3CyBASIS,SyZMAX,SOLN} 
REAL*8 A(M1,N1),0(40),SOLN(20),CBASIS(20),2ZMAX 
INO DAS TO C2ZO0) 5) STEP, S 
COMMON STEP 

M=M1-1 

N=N1-1 

ZERO=1.0D-6 

STEP=0 

DO 1 I=1,M 
CBASIS(I)=C(BASTS(I)) 
CONTINUE 

C(N1)=0.0N0 

DO 2 J=1,N1 

SUM=0.0D0 

DO 3 1=1,™M 
SUM=SUM+ALT,J)*CBASIS(1) 
CONTINUE 
A(M1,J)=SUM-C(J) 
CONTINUE 

DO 4 [=1 aM 
SOLN(I)=AC1I,N1) 
CONTINUE 

ZMAX=A(M1,N1) 

K=0 

XMIN=0.0D0 

DO 8 J=1,N 

X=A(M1,J) 

Px soe ER O)e GOSTO 8 
K=K+1 

IF(X.GE.XMIN) GO TO 8 
XMIN=A(M1,J) 

KIN=J 

CONTINUE 

Eig We Oe OU no 

$=1 

RETURN 

K=0 

XMIN=1.00+20 

DO 10 I=1,M 

X=A(I,KIN) 

PERK tie 20)? GO 10 110 
K=K+1 

X=A(C1sN1)/X 
IF(X.GE.-XMIN) GO TO 10 
XMIN=X 

KOUT=I1 

CONTINUE 

BEC poe Oe TL 

S=2 

RETURN 

K=BASIS(KOUT) 
BASIS(KOUTIJ=KIN 
CBASTS(KOUT)=C(KIN) 
Siew eae eT 

Die (iselctabesiteee 4 a Oia 15 
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RETURN 
15 P=A(KOUT,KIN) | 
DO 13 J=1,Ni1 
A(KOUT »J)=A(KOUT,JI)/P 
13. CONTINUE 
DO 16 {=1,M1 
TPeC USER SKOUT) GO TO716 
P=A(T,KIN) 
DO 14 J=1;,N1 
A(T, J) =A(1,J)-P*#A(KOUT,J) 
14 CONTINUE 
16 CONTINUE 
GO TO 18 
END 
TABLE 


SUBROUTINE TABLETAsN XsY_J5) 
REAL*¥8 Al5yN32)5XsY 
I=N/2 
M=I 

1 FRUXSUEL.AUCU Ms EY) GG TO! S 
TBO SEO. le Oe ome 
I=I1/2 
M=M+]I 
GO Fl Or) 

2 AlLJsM41,2)=ACJI,Mt1,2)+Y 
RETURN 

3) IF(I.NE.~1) GO TO 4 
A(JeMy2)=AlJI_9M,2)4+Y 
RETURN 

4 1l=1/2 
M=M-I 

5 PRUA CS e041 fle et ed) GUT) PT 
PRE Aleteiyi te) TAB? FE) ve; 
I=1/2 
M=M-] 
GO TO 5 ; 

6 AlJoMy2) =A J 2M,2)+Y 
RETURN 

Bf TELIC ED. Woo, Ee 
I=I1/2 
M=M+] 
GO TO l 

8 ACJ sM41,2)=A0I,M+1,.2)4Y 
RETURN 


END 
BOX 


SUBROUTINE BOX(N,R, INT) 

RUA Gert Nia POINT 1234.01 3. 1415926535897 937/,X%/0. 0007 

K=N*¥2 

DO 250 J=lsK 

CALL RANDOM( INT, X) 
aeiay AL ERE IAG @ ko .6 

DO 25 I1=1;N 
22 R(1)=DCOS(2. ODO*PI*POINT (N41) )*(-—2.0D0*DLOG( POINT (1))) **0.5 

RETURN 

END 
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MONTE CARLO METHOD FOR SOLVING 
THE STOCHASTIC LINEAR PROGRAM : 


MAXIMIZE Z, WHERE 


Z= 1.00000 Xl + 2.00000 X2 
SURO Enel oi 
3.200000 Xl + TeOCOOG KA Es 15.00000 
1.00000 Xl + 1.00000 X2 .LE. 10.00090 
AND USING 19000 POINTS 
APPROX. gv ALU Eee yee 1S 20 .02208 
APPROX. VALUE OF VAR(Z) IS 1.03841 
APPROX. VALUE DR Gtx 18) iS Of 0253 
AEPR UX eV ALUGEOReELAZ Jel > 9.99839 
AP PEK sev ALUE, Ole Et AS 15 4.91831 
APPROX. SVALUE OF E(X4) 15 0.00002 


APPROX. VAR-COVARIANCE MATRIX IS 

O25 5140-0 10% 257, 950-0 ls £9 124410800 = —-0.49000D—-056 
eet SOL O52 12550700" =—0 «2020 5D 00 0-23693D-04 
—0.12447D 00 -0.20215D 00 restate) Oe oe ae 
-0.49000D-06 Qe 3o 950-04 §=0.952760 —04 Get le 21-05 
FREQUENCY OF PAIRS MATRIX 


2 3 4 


1 245 0 0 
2 18754 1 
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